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A  Abstract 


Existing  positioning  and  navigation  applications  mainly  rely  on  GNSS.  However,  its  ap¬ 
plication  to  indoor,  metropolitan  and  heavy-foliage  scenarios  is  significantly  constrained 
by  the  satellite  visibility,  limited  accuracy  and  the  intensively  frequency-selective  channel 
condition.  In  this  research,  we  investigated  ranging  and  localization  techniques  using  ultra- 
wideband  (UWB)  signals.  In  particular,  we  have  developed  two  time-of-arrival  (ToA)  esti¬ 
mators  for  single-band  UWB  signals.  To  generalize  them  to  multi-band  UWB  signals,  we 
investigated  the  pros  and  cons  of  coherent  and  non-coherent  multi-band  signal  combining, 
and  studied  the  phase  rotation  calibration  in  order  to  facilitate  effective  coherent  combin¬ 
ing.  Last  but  not  least,  we  developed  cooperative  localization  algorithms  based  on  semi- 
definite  programming  (SDP)  and  developed  algorithms  exploiting  Doppler  shift  in  mobile 
scenarios.  On  these  subjects,  we  have  published /submitted  13  peer-reviewed  journal  and 
conference  papers  [1-13]. 

B  Overview 

The  proposed  research  aims  at  the  fundamental  understanding  and  system-level  devel¬ 
opment  of  alternative  location  and  navigation  in  GPS-denied  scenarios  such  as  indoor, 
metropolitan  and  heavy-foliage  environments,  where  the  application  of  GPS  is  significantly 
constrained  by  the  satellite  visibility,  limited  accuracy  and  the  intensively  frequency-selective 
channel  condition.  When  GPS  is  also  available,  techniques  developed  here  can  be  used  to 
enhance  the  GPS  precision. 

For  time-of-arrival  (ToA)  based  geo-location,  the  precision  of  ToA  estimate  is  crucial.  The 
latter,  however,  is  challenging  in  environments  with  extensive,  and  often  unknown,  multi- 
path,  such  as  heavy  foliage,  indoor  or  inside  caves.  Traditionally,  this  problem  is  treated  as  a 
channel  estimation  problem,  where  the  ToA  is  obtained  as  a  by-product.  However,  such  an 
approach  not  only  induces  significant  complexity  especially  in  extensive  multipath,  but  also 
leads  to  unnecessary  ToA  errors  as  the  optimality  criterion  for  a  channel  estimator  accounts 
for  timing,  amplitude  and  phase  errors  of  all  taps.  To  this  end,  we  developed  estimators  fo¬ 
cusing  on  the  ToA  estimation,  and  without  any  direct  channel  estimation.  Results  show  that 
significant  complexity  reduction  and  accuracy  improvement  can  be  achieved  when  com¬ 
pared  to  the  current  state-of-the-art. 
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For  multi-band  orthogonal  frequency-division  multiplexing  (MB-OFDM)  ultrawi de-band 
(UWB)  systems,  we  estimated  the  ToA  of  wireless  channels  by  using  the  widely  adopted 
equally-spaced  channel  estimator.  By  suppressing  the  energy  leakage  due  to  the  imperfect 
sampling  of  the  leading  channel  paths,  we  proposed  a  novel  ToA  estimator  for  the  MB- 
OFDM  UWB  system.  This  ToA  estimator  is  evaluated  and  compared  with  the  space  alternat¬ 
ing  generalized  expectation  maximization  (SAGE)  algorithm.  Simulation  results  show  that 
the  proposed  ToA  estimator  performs  well  in  all  channel  models  and  outperforms  SAGE  by 
directly  locating  the  first  channel  path. 

We  further  investigated  the  random  phase  rotation  problem  for  MB-OFDM  UWB  sys¬ 
tems.  Random  phase  rotations  need  to  be  calibrated  for  coherent  combining  of  MB  channel 
information  for  high  resolution  time-of-arrival  estimation.  Our  analysis  indicates  that  in  the 
presence  of  random  phase  rotations,  energy  leakage  in  the  estimated  channel  increases  and 
hence  the  ToA  resolution  of  MB-OFDM  signals  degrades.  Based  on  this,  we  propose  to  cal¬ 
ibrate  random  phase  rotations  across  subbands  by  minimizing  energy  leakage.  Simulation 
results  show  that  the  proposed  technique  performs  well  in  802.15.4a  UWB  channels. 

When  it  comes  to  range  measurement  data  fusion  into  positioning  and  location  estimate, 
we  considered  cooperative  localization  by  introducing  and  generalizing  the  application  of 
semidefinite  programming  (SDP).  The  SDP  approach  has  been  widely  used  to  convert  non- 
con  vex  problems  into  convex  ones  in  recent  years.  In  our  work,  we  apply  the  SDP  approach 
to  cooperative  localization  where  the  inter-target  communication  capability  is  exploited  for 
the  purpose  of  coverage  extension  and  accuracy  enhancement.  Cooperative  ToA  and  Re¬ 
ceived  Signal  Strength  (RSS)  minimax  SDP  algorithms  are  also  proposed.  Simulations  show 
that  the  cooperative  localization  with  SDP  can  achieve  satisfying  performance  with  con¬ 
siderably  reduced  complexity.  In  addition,  we  propose  a  virtual  anchor  concept  to  further 
improve  the  localization  accuracy,  especially  in  the  outside-of-the-convex-hull  situations. 

Among  the  various  positioning  techniques,  the  angle-of-arrival  (AoA)  has  traditionally 
been  realized  with  large  and  costly  antenna  arrays.  On  the  other  hand,  motion  of  the  mobile 
target  will  induce  Doppler  which  is  determined  by  the  AoA  of  the  incoming  signal.  Ex¬ 
ploiting  this  Doppler,  we  developed  Doppler-aided  AoA  estimators  for  mobile  targets.  By 
studying  the  geometrical  relationship  among  the  target  and  anchors,  we  successfully  lin¬ 
earized  the  nonlinear  problem,  resulting  in  a  low-complexity  solution.  We  have  also  proved 
the  uniqueness  of  the  solution  using  geometric  tools.  Interestingly,  though  the  problem  is 
formulated  as  AoA-based  geo-location  (since  the  Doppler  can  only  indicate  the  direction  of 
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the  target),  the  solution  is  obtained  via  a  ToA-alike  technique  as  the  position  of  the  target 
is  found  to  be  the  unique  intersection  of  three  circles  (trilateration)  each  determined  by  two 
anchors. 


C  Ranging  with  Single-Band  UWB  Signals:  A  Maximum- 
Likelihood  Approach 

UWB  technology  exhibits  prominent  features  in  many  wireless  communications,  network¬ 
ing  and  localization  applications.  Since  the  ultrashort  pulse  waveform  is  transmitted  at  very 
low  power  in  UWB  systems,  accurate  and  rapid  timing  estimation  becomes  one  of  the  most 
critical  challenges. 

In  this  research,  we  address  this  issue  via  the  establishment  of  a  data-aided  maximum 
likelihood  (ML)  timing  algorithm.  Based  on  the  ML  criterion,  estimation  of  all  multipath 
gains  and  delays  was  pursued  in  some  existing  work  in  the  literature.  However,  they  often 
assume  an  unrealistic  multipath  channel  model  with  no  inter-path  overlapping.  The  real 
channel  with  a  large  number  of  dense  taps  would  make  this  method  impossible  to  imple¬ 
ment.  Here,  we  focus  on  estimation  of  a  single  parameter,  namely,  the  delay  of  the  first 
arriving  path,  without  invoking  any  unrealistic  channel  model  assumption.  We  will  show 
that  our  ML  estimator  is  able  to  collect  multipath  energy,  although  it  does  not  explicitly 
involve  multipath  channel  estimation. 

Considering  the  ML  acquisition  performance  and  the  consistency  requirement,  we  ob¬ 
tain  the  unique  optimum  training  pattern  with  which  the  ML  algorithm  can  be  simplified, 
thus  giving  rise  to  the  simplified  ML  (SML). 

Fine  timing  with  high  accuracy  is  critical  to  localization  with  UWB  technology.  While 
the  data-aided  SML  estimator  can  theoretically  achieve  any  resolution  level,  they  may  suf¬ 
fer  from  the  ambiguity  induced  by  the  weak  tail  of  the  multipath  channel  and  the  extent 
of  the  noise-only  region  between  consecutive  symbols.  To  circumvent  the  ambiguity,  we 
supplement  the  SML  algorithm  with  one  more  step  to  search  for  the  peak  of  the  first-order 
difference  of  the  objective  functions.  By  doing  so,  chip-level  fine  timing  can  be  achieved. 
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C.l  Signal  Model 

In  impulse-radio  UWB  systems,  every  information  symbol  is  transmitted  over  a  duration  of 
Ts  consisting  of  Nj  frames.  During  each  frame  of  Tj  seconds,  a  data-modulated  ultra-short 
pulse  p(t)  with  duration  Tp  «  T/  is  transmitted.  With  binary  pulse  amplitude  modulation 
(PAM),  the  training  symbols  are  drawn  from  the  binary  alphabet  {±1},  and  the  transmitted 
waveform  for  a  single  user  is  modeled  as: 

+OO 

v(t)  =  '/S  Y,  snpT(t  -  nTs ),  (1) 

71=0 

where  £  is  the  energy  per  pulse  and  pr{t)  =  X//=o 1  P(t  ~  jTf  ~  G'Aj  represents  the  symbol- 
long  transmitted  waveform  composed  of  Nj  pulses.  The  pulse  during  the  jth  frame  is  shifted 
by  the  time-hopping  code  cy,  which  takes  integer  value  in  the  range  of  [0,  Nc  —  1],  The  chip 
duration  is  Tc  =  Tf/Nc. 

After  propagating  through  a  multipath  channel  with  Lc  taps,  the  received  waveform  can 
be  written  as  r(t)  =  '  oiiv(t  —  Ti )  +rj(t),  where  a;  and  p  denote  the  attenuation  and  delay 

of  the  /th  channel  tap,  and  77(f)  is  the  zero-mean  additive  white  Gaussian  noise  (AWGN) 
with  power  spectral  density  (PSD)  N0/2.  The  channel  is  assumed  to  be  either  deterministic 
or  quasi-static  over  one  transmission  burst.  We  decouple  the  propagation  delay  r0  from  the 
dispersive  effects  of  the  multipath  channel  by  defining  a  set  of  relative  delays  with  respect 
to  r0,  namely  ti\ 0  =  77  —  r0,  V/.  Without  loss  of  generality  (WLOG),  r0  e  [0,  Ts)  is  assumed  in 
this  work.  Then  the  symbol-long  received  waveform  capturing  the  multipath  channel  effects 
is  given  by: 

Lc~  1 

Pn{t)  =  aiPT<yf  ~  r«|o)  »  (2) 

1=0 

and  the  received  waveform  can  be  rewritten  as: 

+OO 

r(t)  =  VsJ^SnPRjt  -  nTs  -  r0)  +  77(f)  .  (3) 

n=0 

To  develop  the  ML  timing  algorithm,  we  assume  that  ISI  is  absent,  but  inter-frame  interfer¬ 
ence  may  be  present.  This  condition  can  be  easily  satisfied  by  constraining  the  last  frame 
of  each  symbol  such  that  the  nonzero  support  of  pn(t)  does  not  extend  beyond  the  range 
[0,TS).  Note  that  this  setup  can  also  accommodate  high-rate  transmissions  since  the  inter¬ 
frame  interference  is  allowed. 
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For  convenient  manipulation,  we  divide  the  received  signal  into  K  consecutive  symbol- 
long  segments  and  shift  them  so  that  they  all  lie  in  the  range  f  G  [0,  Ts):  rk(t)  =  r(f  + 
/cTs)rect(f),  k  =  1,  •  •  •  ,  K,  where  rect(f)  =  1,  f  e  [0,  Ts),  is  the  window  function.  Substituting 
(3)  and  defining  rjk(t)  =  r/(t  +  kTs)rect(t),  we  have  the  received  segments: 

+oo 

rk(t)  =  Vsy^  snpR(t  +  kTs  -  nTs  —  r0)rect(f)  +  rjk(t),  Vfc  .  (4) 

ra=0 

Since  pR(t)  and  rect(f)  are  confined  within  a  finite  support  [0,  Ts),  it  can  be  easily  induced 
that  for  a  certain  segment  k  only  n  =  k  —  1  and  n  —  k  contribute  nonzero  summands  to  rk(t), 
and  (4)  can  be  explicitly  expressed  as  rk(t)  =  \f£  ( sk-ipR(t  +  Ts  —  r0)  +  skpR(t  —  r0))  rect (t)  + 
pk(t)yk. 

Stack  the  total  K  received  segments  into  a  vector,  and  define  r(i)  =  [r,  (4),  •  •  •  ,  r^(^)]  , 
Si  =  [so,  •  •  •  ,  sK-i]T ,  S2  =  [si,  •  •  •  ,  sk]T  and  17(f)  =  [/71(f),  •  •  •  ,  '/7/r(t)f  •  Then  the  signal  model 
can  be  rewritten  in  the  following  compact  vector  form: 

r(t)  =  VSsip^\t-,  t0)  +  VSs2pR\t ;  r0)  +  r/(f)  ,  (5) 

where  p^(f;r0)  =  p_R(f  +  Ts  —  r0)rect(f)  and  p^it:  r0)  =  pR(f  —  r0)rect(f)  consist  of  the 
circularly  shifted  version  of  the  symbol-long  received  waveform  pR(t).  It  is  noteworthy  that 
PR\t]  To)  and  pR\t;  r0)  are  not  overlapping  in  time;  that  is,  the  former  is  strictly  zero  for 
f  6  [r0,  T,),  and  the  latter  is  strictly  zero  for  f  e  [0,  To). 

C.2  ML  Timing  Algorithm  and  Its  Acquisition  Performance 

In  this  section,  we  will  first  develop  the  ML  timing  algorithm  for  arbitrary  known  transmit¬ 
ted  symbol  sequences,  and  then  evaluate  the  timing  acquisition  performance  of  the  algo¬ 
rithm. 

C.2.1  ML  Timing  Algorithm 

With  the  signal  model  in  (5),  the  deterministic  but  unknown  parameters  are:  i)  the  overall  re¬ 
ceived  symbol-long  waveform  pR{t)  (or  equivalently,  its  circularly  shifted  version  pR\t]  r0) 
and  pR\t ;  t0))  which  carries  the  dispersive  multipath  channel  information;  and  ii)  the  prop¬ 
agation  delay  t0.  Given  pR(t)  and  r0,  the  log-likelihood  function  for  (5)  bears  the  form: 

rTs 

In  A  (r(f);pfl(f),  To)oc  /  -||  r(f)  -  V^Sip^t;  t0)  -  V£s2p^\t]  r0)  ||2  dt 
Jo 
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c xj  2 VSrT(t)  (^siPft\t]  T0)  +  s2p(*\t ;  r0))  -5  ||  r0)  +  s2pj)(t;  r0)  ||2dt.  (6) 

Our  task  is  to  obtain  the  ML  estimates  for  p#(f)  and  r0  by  maximizing  (6).  We  use  the 
notation  x  to  indicate  a  conjecture  of  unknown  parameter  x.  The  ML  estimation  will  be 
accomplished  in  two  stages:  based  on  a  fixed  conjecture  t0/  we  first  obtain  pR(t ;  r0)  as  a 
function  of  r0;  then  we  replace  pR(t)  with  pR(t:  t0)  in  (6)  to  find  the  ML  estimate  of  f0. 

In  the  first  stage,  the  integral  can  be  removed  without  affecting  the  optimality,  since  the 
ML  estimate  of  pR(t)  is  to  be  obtained  in  an  instantaneous  manner.  As  emphasized  before, 
given  a  candidate  f0,  pR\t]  f0)  and  pR\t ;  f0)  are  non-overlapping  in  time.  Accordingly,  we 
divide  r(t)  into  two  disjoint  parts  in  time  as  well:  r(f)rect(f  +  Ts  —  f0)  for  t  G  [0,  f0)  and 
r(f)rect(f  —  f0)  for  £  G  [f0,  Ts).  The  circularly  shifted  waveforms  p^(£;  f0)  and  pR\t;  f0)  can 
thus  be  estimated  separately.  Specifically,  the  objective  function  for  pR\t',  r0)  is  [c.f.  (6)]: 

J(a)(f;f0)  =  2v/£rT(f)sip^)(t;f0)  -£  ||  Sip^fjfo)  ||2,  t  G  [0,to)  . 

Taking  the  derivative  of  J^(f;  f0)  with  respect  to  the  instantaneous  pR\t;  f0),  and  setting  it 
to  zero,  we  have  the  ML  estimate  of  (£;  f0): 

1  A' 

Pfl  (i;  ^o)  =  — y=  Sk-irk(t)rect(t  +  TS-  t0),  £  G  [0,  r0)  .  (7) 

v  fc=i 

Likewise,  the  ML  estimate  of  pR\t;  tq)  can  be  obtained  by  maximizing  the  following  objec¬ 
tive  function  [c.f.  (6)]: 

J(h\t]f0)  =  2VSrT(t)s2p^\t]f0)  -S  ||  s2p^\t-,f0)  ||2,  £G  [t0, Ts)  , 

and  the  resultant  estimate  is: 

1  K 

PR\t] To)  =  XJ  skrk(t)rect(t  -  t0),  t  G  [t0,  Ts)  .  (8) 

In  the  second  stage,  we  plug  (7)  and  (8)  back  into  (6),  and  discard  the  norm  square  term 
whose  integral  is  not  affected  by  the  candidate  f0.  It  turns  out  that  the  new  ML  objective 
function  for  r0  becomes: 

K  I<  \ 

rT(t)si  ^  sk-irk(t)rect(t  +  TS-  f0)  +  rT{t)s2  ^  skrk{t)rect{t  -  r0)  1  dt 

k= 1  k= 1  / 

X  K  K  rTs 

=-pEE/„  (t fc(f)'Sm-iSfc_irect(f+Ts  — r0)+r m{t)r k(t)smskrect(t  —  T0))  dt,  (9) 

m=  1  k=  1'  0 


jUL{fo)  =  ^l 
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and  the  ML  estimate  of  r0  can  be  obtained  by  maximizing  JML(ro): 

f0  =  argmax  JML(f0)  .  (10) 

TO 

Proposition  1  [ML  Timing  Estimation]  The  ML  timing  estimator  can  be  implemented  in  four 
steps: 

•  Step  1:  Take  K  received  segments  rk(t),  k  =  1, 2,  •  •  •  ,  K; 

•  Step  2:  Tor  each  candidate  f0  calcidate  K 2  cross  (and  auto)  correlations  among  all  pairs  of  the 
segments  as  in  (9); 

•  Step  3:  Average  the  K2  correlations  as  suggested  by  (9); 

•  Step  4:  Choose  the  f0  which  maximizes  Jml(t0)  as  the  ML  estimate  f0  according  to  (10). 

From  Proposition  1  one  should  be  aware  that  the  computational  complexity  of  the  ML 
timing  estimator  is  very  high.  For  each  fo  candidate,  one  need  to  calculate  K2  correlations 
and  K 2  summations.  The  high  complexity  is  expected  to  be  reduced  for  practical  imple¬ 
mentation. 


C.2.2  Acquisition  Performance  of  the  ML  timing  Algorithm 


In  order  to  evaluate  the  performance  of  the  ML  estimator,  we  need  to  find  the  statistical 
properties  of  the  objective  function.  Re-express  JML(f0)  as  the  sum  of  its  noise-free  part  and 
noise  term  as  JML(f0)  =  +  £ml(t0).  Let  us  first  consider  the  noise-free  part 

1  K  [Ts 

{Pm(t)Pk(t)sm-isk-irect(t  +  Ts  -  T0)+pm(t)pk(t)smskrect(t  -  t0))  dt  (11) 

m,k=T° 


where  pk(t)  denotes  the  signal  part  of  rk(t).  Assuming  the  candidate  r0  <  r0  WLOG, 

we  further  define  Sai  —  VS  JqS~T0 p2R(t)dt,  £a2(t0)  =  y/S  fR3ff°+w  p2R(t)dt,  and  £b(tq)  — 
'/£  ffss_To+foP2R{t)dt.  Then  (11)  becomes 


</oML(ro) 


1 

lc2 


K 

(SmSk^Al  +  Sm-lSk-l^A2(To)  +  Srn_iSmSfe_iSfc£s(fo)) 

m,k= 1 


l 

K2 


K 

E  (^(^°)  + 

m,k= 1 


Sm—  lSmSk-lSk£B(fo)) 


(12) 
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where  £a(t0)  =  £ai+£a2(tq)-  Notice  that  £a{tq)+£b{tq)  =  £  Jo*  pji(t)dt  =  Sr  is  the  unknown 
but  constant  energy  of  a  received  segment  independent  of  the  trial  value  fo-  The  noise-free 
part  of  the  objective  function  can  thus  be  simplified  as 

J0ML(f0)  =  £r  -  §ffS(f0)  ,  (13) 

where  AK  =  Y^m,k= i(l  -  Sm-ismSk-iSk )  is  a  positive  (as  long  as  not  all 
1,  Vm,  k,  which  can  be  easily  avoided.)  parameter  determined  by  the  transmitted  pilot  se¬ 
quence.  By  definition,  the  condition  of  correct  timing  f0  =  r0  ensures  that  SB( f0)  vanishes 
and  Jo1L(fo)  achieves  its  unique  maximum  Sr.  Notice  that  since  A K  is  a  sum  of  K 2  constants, 
its  value  is  on  the  order  of  K2,  or  can  be  explicitly  written  as  c\  K2,  where  a  is  a  constant. 

We  now  go  to  the  noise  term  £ML(f0).  It  can  be  proved  that  £ML(f0)  is  Gaussian  distributed 
with  mean  and  variance: 


mMLm  =  ^ 


var{£ML(f0)}  = 


2-/Vp  Jq  (fo) 
K 


+ 


n2bts 

K 2 


(14) 


Then,  one  can  obtain  the  mean  and  variance  of  the  overall  Gaussian  ML  objective  function: 


E{JML(f0)} 
var{  JML(f0)} 


+  E{«ML(f„)}  =  ./0ML(f „)  +  ^ 
var™}  =  2^|A)+^. 


(15) 


As  one  can  see  from  (15),  JML(f0)  asymptotically  converges  to  ./()1L(f0)  as  K  -a-  oo,  suggest¬ 
ing  the  optimality  of  the  ML  estimator. 

We  adopt  the  probability  of  detection  lower  bound  Pd  to  evaluate  the  coarse  timing 
(acquisition)  performance  of  the  ML  algorithm.  Instead  of  estimating  the  true  r0,  coarse 
timing  aims  at  finding  no  such  that  n0Ti  —  r0|  <  Tu  where  Tr  is  the  searching  step  size  in  the 
ML  algorithm.  Correspondingly,  the  maximization  problem  in  (10)  becomes: 


h0  —  argmax  JML(h0Ti)  (16) 

no 

and  the  probability  of  detection  is  given  by 


pML  =  Y>r{ho  =  n0}  =  Prjmax  JML(h0T,:)  =  JML(n0Ti)}  .  (17) 

no 

Since  JML(n0Ti )  is  Gaussian  distributed,  the  lower  bound  of  (17)  is  given  by 

pML  _  T-r  F  (  E{JML(noT,)}-E{JML(h0rt)}  \ 
no}no  V  Vvar{jML(no Ti)}  +  var{ JML(hoTi)}  ) 
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where  F(-)  is  the  cumulative  distribution  function  (cdf)  of  Gaussian  distribution  with  zero 
mean  and  unit  variance. 

Substituting  the  mean  and  variance  of  the  objective  function  given  in  (15),  we  can  obtain 
the  probability  of  detection  lower  bound: 


pML 
±-  d 


AA'^(noTj) 


s/lNQ£RK*-2NQ&K£B{h0Ti)K+2N%BTaKH 


(19) 


Remarks:  i)  As  K  increases  the  variance  of  the  objective  function  decreases  and  the  proba¬ 
bility  of  detection  lower  bound  increases.  This  suggests  that  the  timing  performance  would 
benefit  from  more  correlation  averaging;  and  ii)  As  A k  increases  the  variance  of  the  objec¬ 
tive  function  is  reduced  and  the  probability  of  detection  lower  bound  increases.  Intuitively, 
(13)  provides  another  evidence  that  the  objective  function  becomes  sharper  along  with  the 
increase  of  A  k,  making  detection  of  no  easier.  Since  A  k  is  determined  by  the  transmitted 
symbol  sequence,  one  can  expect  that  the  acquisition  performance  would  be  markedly  im¬ 
proved  by  optimizing  the  training  sequence  such  that  A K  is  maximized. 


C.3  Training  Sequence  Design  and  The  SML  Algorithm 

C.3.1  The  Optimum  Training  Sequence  Pattern 

The  value  of  AK  =  fc=i(l  —  sm-ismSk-iSk)  is  determined  by  the  signs  of  the  consecutive 
symbols.  Define  Ck  —  Sk-iSk  E  (±l},/c  =  l,---  , K,  as  the  product  of  two  consecutive 
symbols.  Then  Ck  belongs  to  one  of  the  two  groups,  T+  A  [ck.  VA;  :  Ck  =  1}  with  cardinality 
K. |_  and  T_  A  {cfc,  V/c  :  Ck  =  —1}  with  cardinality  K_.  Evidently,  K+  +  K_  =  K. 

Lemma  1  For  a  specific  I\ ,  max{  AK}  =  K2  is  achieved  when  K+  =  K_  =  K/ 2;  that  is,  half  of 
the  {ck}k  elements  belong  to  T+  and  the  other  half  belong  to  T_. 

Notice  that  we  only  considered  even  K  WLOG,  since  odd  K  has  the  same  maximization 
result  with  its  even  neighbor  K  —  1. 

Lemma  1  gives  the  condition  that  maximizes  A k  for  a  particular  K .  Additionally,  the 
optimum  training  sequence  should  also  be  consistent;  that  is,  applicable  to  arbitrary  K .  To 
this  end,  we  first  notice  that  the  ML  timing  estimator  requires  K  >  2,  since  when  K  =  1, 
J0ML(f0)  =  £r  [c.f.  (12)]  is  simply  a  constant  and  gives  no  information  about  r0.  hollowing 
Lemma  1,  the  consistency  property  can  be  ensured  VA'  >  2  by:  i)  partitioning  the  {ck}k  se¬ 
quence  into  doublets  {c2n-i  c2n},  n—  1,  2,  •  •  •  ,  K/2;  and  ii)  making  sure  that  each  {c2n_i  c2n} 
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doublet  contains  a  “+1"  and  a  “—1",  i.e.,each  doublet  should  be  either  {+1,  —1}  or  {—1,  +1}. 
With  this  condition,  the  {ck}k  sequence  always  has  K+  =  K_  regardless  of  K.  Note  that,  the 
K  =  2  case  which  ensures  rapid  acquisition  using  only  2  segments  is  a  natural  corollary  of 
the  consistency  property 

In  addition,  to  guarantee  that  the  “+1” ,  "—1"  pairing  condition  holds  for  any  doublet 
starting  from  odd-  and  even-indexed  symbols,  all  the  doublets  should  be  identical.  In  other 
words,  they  are  either  all  {+1,  —1}  or  all  {—1,  +1}.  As  a  result,  this  gives  rise  to  a  unique 
training  sequence  {sk}k  consisting  of  the  repeated  pattern  {+1,  +1,  —l,  —1}  (or  its  circularly 
shifted  versions).  We  summarize  the  analysis  in  the  following  result: 

Proposition  2  [Optimum  Training  Sequence]  The  unique  optimum  training  sequence  for  the 
ML  estimator  has  the  structure 

sk  =  (~ l)Lfc/2j  ,  (20) 

which  ensures  rapid  acquisition  using  as  few  as  3  symbols  and  is  applicable  to  arbitrary  K(>  2). 

C.3.2  Simplified  ML  (SML)  Algorithm 

For  simplicity,  denote  the  integrals  in  the  objective  function  for  r0  (9)  as  Considering 
the  partition  of  the  training  sequence  by  groups  V  ,  and  T_,  we  can  rewrite  (9)  as 

dML(fo)  =  -^2  <  .W  +  .W  (21) 

^  {(m,k):Cm,Ck&T+}  {(?n,fc):cm,cfce T_} 

+  ^  '  jm,k  +  ^  ^  jm,k 

{(m,fc):cmer+,cfcer_}  {(m,fe):cmer_,cfcer+} 

Consider  the  first  two  summations.  Since  cm  and  cn  are  chosen  from  the  same  group  (namely 
T+  in  the  first  summation  and  T_  in  the  second  summation),  the  noise-free  parts  of  the 
summands  are  exclusively  £a(to)  +  £b(t0)  =  Sr  [c.f.(12)].  Furthermore,  it  is  not  difficult 
to  verify  that  the  noise  terms  in  the  first  two  summations  do  not  change  with  the  shift 
candidate  f0.  Therefore,  the  first  two  summations  are  nothing  but  constants,  which  provide 
no  information  on  r0.  If  one  knows  which  (cm,Ck)  pairs  give  rise  to  these  summands,  one 
can  avoid  calculating  their  corresponding  cross  correlations. 

The  optimum  training  sequence  given  by  (20)  precisely  allows  one  to  achieve  this.  The 
repeated  pattern  {+1,+1,—  1,  —  1}  indicates  that  the  received  K  symbol-long  segments  can 
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be  divided  into  two  groups  by  simply  checking  their  indices.  Specifically,  if  the  symbol- 
long  received  segment  with  odd  index  r2k_i(t)  carries  two  successive  symbols  satisfying 
c2k- 1  =  S2k-2$2k-i  €  r+  (or  T_),  then  the  symbol-long  received  segment  with  even  index 
r2k  must  carry  two  successive  symbols  satisfying  c2k  =  s2k_±s2k  G  T_  (or  T+).  Retaining 
only  the  cross  correlations  between  the  two  groups,  we  can  obtain  the  simplified  ML  (SML) 
objective  function  as 

9  K/2 


Jo 

+r2m-it)r2k{t)  (-1)  rect (t  -  t0)  dt .  (22) 

Note  that  the  last  two  summation  terms  in  (21)  are  exactly  the  same,  which  explains  the 
reason  why  a  scaling  coefficient  2  shows  up  in  (22). 

We  can  rewrite  JSML(f0)  by  putting  the  double  summations  into  the  integral: 

1  rTs  (  2  K/2  \  (  2  K'2 

=  2/  Ly(-i)%-,w  Lyi-itSiW 

y  m=  1  /  \  V  k=1 

■  (rect(f  +  TS  —  f0)  —  rect(f  —  f0))  dt  .  (23) 

The  above  integrand  includes  the  product  of  three  terms.  The  first  is  the  average  of  the  odd 
indexed  received  segments  which  satisfies  the  condition  for  the  group  T+;  the  second  is  the 
average  of  those  even  indexed  received  segments  falling  into  the  group  T_;  and  the  last 
term  is  the  window  function  accounting  for  the  guess  shift  f0. 

Proposition  3  [SML  Timing  Estimation]  By  employing  the  optimum  training  sequence  given 
in  Proposition  2,  the  SML  estimator  can  be  implemented  with  much  lower  complexity  than  the  ML 
estimator: 

•  Step  1:  Take  K  received  segments  rk(t),\/k; 

•  Step  2:  Average  the  odd  and  even  indexed  segments  respectively  as  suggested  by  (23); 

•  Step  3:  For  each  candidate  f0,form  the  window  functions  rect(t  +  Ts  —  f0)  and  —rect(t  —  f0), 
and  calcidate  Jsml(t0 )  as  (23); 
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•  Step  4:  Choose  the  r0  which  maximizes  Jsml(t0 )  as  the  SML  estimate  f0;  that  is,  f0  = 
argmaxf0  JSML{jo). 


It  is  important  to  note  that  the  complexity  of  the  SML  timing  estimator  is  significantly 
reduced.  One  only  needs  to  evaluate  K  summations  and  1  correlation  for  each  r0  candi¬ 
date.  Moreover,  in  a  digital  implementation,  we  do  not  need  to  compute  the  correlation  for 
every  new  f0,  as  most  of  the  correlation  is  identical  from  the  current  f0  value  to  the  next. 
Additional  computing  saving  can  be  obtained  by  only  updating  the  difference  instead  of 
calculating  every  correlation  anew. 

Like  the  ML  estimator,  we  are  also  interested  in  the  acquisition  performance  of  the  SML 
estimator.  Inherited  from  (12)  with  s2m-2S2m-iS2k-iS2k  =  — 1 ,  Vm,  k,  the  noise-free  part  of  the 
SML  objective  function  in  (22)  can  be  expressed  as: 


jSML(~^  _  £4(70)  £b(tq) 


Sr  —  2£b(to) 
2 


(24) 


The  noise  term  £SML(r0)  is  Gaussian  distributed  with  zero  mean  and  variance  (the  proof  is 
similar  to  that  for  the  ML,  thus  omitted  here): 


var{£SML(f0)} 


N0Sr  NyBTs 
2  K  +  2K2 


(25) 


After  calculating  the  mean  and  variance  for  JSML(r0): 

E{ JSML(f0)}  =  £Ah)  =  £fl(fo)  ,  var{ JSML(fo)} 


Nq£r  NlBT, 

2  K  2K2  ’ 


(26) 


we  obtain  the  probability  of  detection  lower  bound  for  the  SML  algorithm  as: 


pSML 

f-d 


no^no 


K£B{n0Ti ) 


\/  KNqSr  +  N2BTS 


(27) 


C.4  Implementation  Considerations  and  Simulated  Performance 

Theoretically,  the  SML  algorithm  can  always  detect  n0  such  that  \n{)Tt  —  r0|  <  T,  on  any 
resolution  level  X \  as  long  as  the  complexity  of  the  receiver  is  allowed.  In  practice,  however, 
significant  attenuation  at  the  tail  of  a  multipath  channel  and  the  extent  of  noise-only  region 
between  consecutive  symbols  make  things  more  complicated.  The  unique  peak  at  n0  of  the 
objective  function  tends  to  be  comparable  with  its  left  neighbors  h0  <  "0  even  when  the 
noise  is  absent.  Therefore,  the  correct  timing  of  ho  =  no  with  fine  (e.g.,  chip-level)  resolution 
is  not  easily  distinguishable  as  the  peak.  On  the  other  hand,  we  notice  that  the  value  of  the 
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objective  function  decreases  dramatically  for  n0  >  n0r  as  a  result  of  the  first  few  strong  taps 
of  the  channel.  Thanks  to  the  different  behavior  of  the  regions  n0  <  no  and  n0  >  no,  one  can 
resort  to  the  first-order  difference  of  the  objective  function  in  aid  of  finding  n0  at  chip  level. 

Suppose  that  the  frame-level  acquisition  has  already  been  achieved.  After  obtaining 
the  values  of  J(h0Tc )  in  the  right  frame  where  n0  is  located,  take  the  difference  A  J(n0Tc)  = 
J(h0Tc)  —  J(h0Tc  —  wTc),  where  w  G  [1,  Nc]  and  wTc  denotes  the  step  size.  Then  the  candidate 
ho  which  maximizes  A  J(h0Tc)  will  be  regarded  as  the  estimate  of  n0. 

We  use  the  channel  model  IEEE  802.15.3a  CM1  to  generate  the  multipath  channel.  The 
UWB  pulse  is  the  second  derivative  of  the  Gaussian  function  with  unit  energy  and  duration 
Tp  «  1  ns.  The  frame  duration  is  Tf  =  35  ns,  and  each  symbol  contains  Nf  =  32  frames.  A 
random  time  hopping  code  Cj  is  uniformly  distributed  over  [0,  Nc  —  1],  with  Nc  =  35  and 
Tc  —  1  ns.  To  avoid  ISI,  the  time  hopping  code  for  the  last  frame  of  each  symbol  is  set  to 
CNf- 1  =  0. 

Test  1  illustrates  the  optimality  of  the  training  sequence  given  in  (20)  for  the  ML  estima¬ 
tor.  Fig.  1(a)  compares  the  optimum  training  pattern  with  four  randomly  chosen  sequences. 
We  can  observe  that  for  K  =  2  only  the  optimum  training  sequence  works.  With  K  >  2,  al¬ 
though  the  other  four  can  also  work,  the  optimum  training  sequence  consistently  provides 
the  best  performance  for  any  K.  We  will  employ  the  optimum  training  sequence  in  the 
subsequent  simulations. 

Test  2  plots  the  objective  functions  of  the  ML  and  SML  algorithms  in  one  realization  at 
the  frame  level  in  Fig.  1(b).  For  comparison,  the  corresponding  noise-free  part  of  the  SML 
(and  TDT)  algorithm  is  also  provided.  It  is  shown  that  the  ML  and  SML  objective  functions 
have  identical  shape.  The  difference  between  them  remains  the  same  for  all  candidates  ho- 

Test  3  depicts  the  acquisition  performance  of  the  SML  algorithm.  The  probabilities  of 
detection,  together  with  the  analytical  lower  bounds  (27)  are  plotted  in  Fig.  2(a).  The  nor¬ 
malized  mean  square  errors  (MSE),  which  are  normalized  with  respect  to  T2,  are  also  shown 
in  Fig.  2(b). 

Test  4  is  designed  to  show  the  chip-level  fine  timing  performance.  The  frame-level  acqui¬ 
sition  is  assumed  to  be  achieved  beforehand.  As  illustrated  in  Fig.  3(a),  the  performance  of 
the  difference  operation  depends  on  the  step  size  wTc.  Accordingly,  we  choose  the  optimum 
value  for  CM1  in  our  simulation  as  wTc  =  3ns  at  low  SNR  and  wTc  =  8ns  at  high  SNR.  The 
normalized  MSE  for  the  SML  algorithm  with  various  K  is  plotted  in  Fig.  3(b).  Notice  that 
all  curves  reach  an  error  floor  since  the  timing  with  chip-level  resolution  is  performed. 


13 


(a) 


(b) 


Figure  1:  (a)Pd  for  various  sequences,  £/Nq  =  ldB.  (b)  Objective  function  magnitude  for  the  ML 
and  SML  algorithms,  £/Nq  =  5 dB,  K  =  8  and  no  =  15. 


(a) 


(b) 


Figure  2:  (a)  Pd  and  Pd  for  the  SML  algorithm,  coarse  timing,  (b)  Normalized  MSE  for  the  SML 
algorithm,  coarse  timing. 
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(a)  (b) 

Figure  3:  (a)  Normalized  MSE  versus  differential  step  wTc  under  various  £/Nq.  (b)  Normalized 
MSE  for  the  SML  algorithm,  fine  timing. 

C.5  Summary 

In  this  research,  we  developed  the  data-aided  ML  timing  algorithm,  and  derived  the  opti¬ 
mum  training  sequence  for  the  ML  algorithm.  Based  on  this  optimum  sequence,  the  original 
ML  algorithm  can  be  simplified  without  affecting  its  optimality.  Extensive  simulations  have 
been  performed  to  corroborate  our  theoretical  analysis. 

D  Ranging  with  Single-Band  UWB  Signals:  Focusing-on- 
First-arrival  (FoFa) 

In  the  literature,  various  techniques  have  been  proposed  for  ToA  estimation  for  broadband 
wireless  systems.  A  straightforward  method  is  to  capture  the  first  signal  arrival  by  a  slid¬ 
ing  correlator  with  a  locally  generated  template  at  the  receiver.  In  order  to  achieve  high- 
accuracy  channel  estimation  and  consequently  high-accuracy  ToA  estimation,  the  correlator 
has  to  slide  at  a  sufficiently  small  step  size.  This  implies  either  very  long  search  time  or  very 
high  complexity  receiver,  depending  on  whether  the  correlations  are  performed  in  serial  or 
parallel. 

In  order  to  avoid  the  sliding  correlator  required  in  the  time  domain  approaches,  fre¬ 
quency  domain  alternatives  were  proposed  for  ToA  estimation,  which  are  typically  based 
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on  the  frequency  domain  equalization  (FDE)  framework.  This  type  of  techniques  treat  ToA 
estimation  as  a  by-product  of  channel  estimation.  Intuitively,  if  one  obtains  all  channel  pa¬ 
rameters  including  each  path's  gain  and  delay,  then  the  delay  of  the  first  multipath  arrival 
automatically  becomes  the  ToA  estimate.  Using  the  fast  Fourier  transform  (FFT),  FDE  can 
be  realized  at  a  low  complexity.  For  FDE  systems,  the  frequency  domain  channel  informa¬ 
tion  should  be  converted  to  the  time  domain  to  facilitate  ToA  estimation.  This  conversion 
gives  rise  to  various  channel  estimation  techniques  including  the  model-based  estimator 
that  models  the  channel  as  a  tapped  delay  line  structure,  the  space  alternating  generalized 
expectation  maximization  (SAGE)  algorithm  and  the  subspace-based  estimator. 

Problems  arise  when  these  channel  estimators  are  employed  for  ToA  estimation.  Firstly, 
the  fixed  structure  of  the  tapped  delay  line  model  can  not  match  the  contiguously  changing 
tap  delays.  Secondly,  the  criterion  used  for  channel  estimation  purpose  is  often  either  to 
maximize  the  likelihood  or  to  minimize  the  estimation  error  between  the  estimates  of  all 
channel  paths  and  the  real  ones.  As  a  result,  the  contribution  from  the  first  path  is  very 
small,  especially  when  there  is  a  large  number  of  dense  multipath  components.  Hence, 
these  approaches  often  encounter  poor  convergence  or  suffer  from  local  optimum  problems, 
when  being  used  as  ToA  estimators. 

To  avoid  these  problems,  we  propose  to  estimate  the  ToA  of  the  channel  by  Focusing  On 
the  First  Arrival,  which  we  term  as  FoFa.  Compared  to  conventional  methods,  the  FoFa  ToA 
estimator  has  the  following  advantages.  Firstly,  FoFa  uses  the  frequency  domain  channel 
information  readily  available  in  FDE  systems  to  recover  the  time  domain  channel  with  an 
equally-spaced  model,  which  can  effectively  reduce  the  complexity  of  the  estimation.  Sec¬ 
ondly,  unlike  conventional  ToA  estimators,  FoFa  relies  on  a  novel  criterion  to  optimize  the 
ToA  estimation;  that  is,  locating  the  first  path  by  minimizing  the  energy  leakage  prior  to  the 
first  path.  As  a  result,  FoFa  directly  addresses  the  delay  estimation  of  the  first  channel  path, 
without  unnecessarily  caring  about  the  estimation  errors  for  the  trailing  paths.  This  avoids 
the  poor  convergence  and  local  optimum  problems  often  encountered  when  treating  ToA 
estimate  as  a  by-product  of  channel  estimation.  In  addition,  all  computations  are  carried 
out  with  the  baseband  signal  instead  of  the  analog  or  oversampled  waveform.  Therefore, 
the  complexity  of  the  FoFa  ToA  estimator  is  much  lower  than  that  of  the  sliding-correlator- 
based  estimator. 

Notation:  We  will  use  bold  upper  and  lower  cases  to  denote  matrices  and  column  vectors, 
respectively.  We  will  use  (-)r  and  (-)n  for  transpose  and  conjugate  transpose  of  matrices 
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Figure  4:  Baseband  transceiver  diagrams  for  OFDM  and  SC-FDE. 
and  vectors,  and  (■)*  for  conjugate  of  complex  numbers. 

D.l  Broadband  System  Model  with  FDE 

In  this  section,  we  will  briefly  review  the  FDE  systems.  Typical  FDE  systems  include  the 
well-known  orthogonal  frequency-division  multiplexing  (OFDM)  system  and  the  single¬ 
carrier  frequency  domain  equalization  (SC-FDE)  system.  Let  us  first  explain  the  basic  idea 
of  OFDM  systems,  and  then  briefly  introduce  SC-FDE.  At  the  OFDM  transmitter  (see  Fig. 
4),  a  block  of  information  symbols  s  =  [si, . . . ,  sn]t  are  multi-carrier  modulated  onto  N  or¬ 
thogonal  digital  subcarriers  to  form  x  =  FHs,  where  F  is  the  FFT  matrix.  A  guard  interval 
(GI)  in  the  form  of  padding  zeros  (ZP)  or  cyclic  prefix  (CP)  is  added  to  each  block  to  mitigate 
the  inter-symbol  interference  (ISI).  After  the  digital-to-analog  conversion  (DAC),  the  signal 
is  carrier  modulated  and  transmitted  from  the  antenna.  The  transmitted  signal  then  propa¬ 
gate  through  the  channel:  h(t)  =  Ytl= i  hr  5(t  —  77),  where  {hi}f=1  and  {t'Z}^=1  are  amplitudes 
and  delays  of  the  L  channel  paths,  respectively.  Note  that  we  do  not  require  {r/}(=l  to  be 
equally  spaced  and  do  not  assume  any  a  priori  channel  information. 

At  the  receiver,  the  arriving  waveform  is  carrier  demodulated  and  analog-to-digital 
(A/D)  converted  to  baseband  discrete-time  samples.  After  the  GI  is  removed,  the  base¬ 
band  signal  is  multi-carrier  demodulated  with  FFT  operation  to  generate  sequence  {rk}k=1- 
With  some  coarse  synchronization,  it  can  be  easily  shown  that 

rk  =  skHk  +  wk  ,  k  =  1, . . . ,  N  ,  (28) 

where  wk  is  the  additive  white  Gaussian  noise  (AWGN),  and  {Hk}k=1  are  the  Fourier  trans¬ 
form  (FT)  coefficients  of  the  channel: 

L 

Hk  =  ^hr  exp {-jukn)  ,  uk  =  2kn/TB  ,  (29) 

i=i 
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where  Tb  is  the  information  block  duration.  Based  on  (28),  the  maximum  likelihood  (ML) 
estimate  of  Hk  can  be  formed  as  Hk  =  rkjsk,  1  <  k  <  N. 

Fig.  4  compares  the  baseband  block  diagrams  of  OFDM  and  SC-FDE.  The  difference 
between  OFDM  and  SC-FDE  is  where  the  IFFT  operation  is  performed.  In  OFDM,  IFFT 
is  adopted  at  the  transmitter  to  modulate  information  symbols  on  subcarriers.  In  SC-FDE, 
IFFT  is  performed  at  the  receiver  to  convert  the  FDE  output  back  to  the  time  domain  sym¬ 
bols. 

D.2  Energy  Leakage  in  Channel  Estimate 

Following  the  literature,  the  channel  is  estimated  by  fitting  the  equally-spaced  model  (30) 
to  the  frequency  domain  channel 

L 

h(t)  =  y ^Jin-8(t-  fn)  (30) 

71=1 

where  fn  =  f\  +  (n  —  l)Tp/  1  <  n  <  L  and  L  =  [~T/,/T/;]  with  Th  being  the  maximum 
channel  delay  spread  and  Tp  the  tap  interval  which  is  chosen  as  the  inverse  of  the  signal 
bandwidth.  By  doing  this,  ToA  of  the  channel  can  be  estimated  by  optimizing  the  single 
free  parameter  f\.  This  is  in  comparison  with  existing  channel  estimators  which  search  all 
individual  channel  path  delays  in  [0,  Th\. 

Given  flr  {hn(fi)}^=1  are  expected  to  satisfy  [c.f.  (29)] 

L 

^  ^  ^7i(Tl)6X p(  J0Jb,k'^n)  k  (31) 

71=1 

where  the  dependance  of  {hn{fi)}n=i  on  A  is  explicitly  shown.  More  compactly  written,  we 
have 


G{f\)h(fi)  =  H  +  ri  (32) 

where  H  =  \H±,  Hi , . . . ,  HTB~\T  and  ij  =  [rjJ,  . . . ,  rj^]1  include  all  the  subband  FT 
coefficients  and  noise  terms,  respectively,  h{ji)  =  ^2(A),  •  •  • ,  hi(fi)]r  and  G(fi)  = 

\G\  (fi),  Gj2  (fi), . . . ,  Gb(ti)]T  with  Gb(fi)  being  a  K  xL  FT  matrix  with  the  (k,  n)th  element 
being  exp(—  jujbjkfn).  Based  on  Eq.  (32),  h(f\)  can  be  obtained  by 

Hn)  =  (Gn(f1)G(f1))-1  Gh(t1)H  +  ft  (33) 
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Figure  5:  (a)  Strong  energy  leakage,  (b)  Weak  energy  leakage. 


where  fj  is  the  noise. 

From  (29),  the  frequency  domain  channel  response  Hbtk  contains  information  from  all  L 
channel  paths.  Therefore,  the  channel  estimate  can  also  be  expressed  as 

L 

h(f i)  =  Y.  Gn(f1)H(l)  +  r,  (34) 

1=1 

where  H{1 )  is  the  frequency  domain  contribution  from  the  Zth  channel  path  which  can  con¬ 
structed  by  H(l)  =  [H^(l),Ho{l),.  ■  ■  ,H1j3(l)]T  andHb(l)  =  [hiexp(-jub)1Ti),  hiexp(-jubt2Ti), . . 
hiexp(—juJb,KTi)]T.  Same  as  other  ToA  estimators  for  MB-OFDM,  Eq.  (33)  also  assumes  that 
no  random  phase  rotation  exists  in  subband  signals. 

Proposition  4  For  the  Ith  channel  path  ( hi ,  77),  given  the  channel  estimator  (33),  only  the  mth  tap  of 
the  channel  estimate  contains  non-zero  contribution  from  the  Ith  path  if  this  path  is  exactly  sampled 
by  the  mth  tap  as  ti  =  rm,  3m  e  [1,  L\. 

This  is  because  when  77  =  fm/  H (7)  will  be  the  mth  column  of  the  matrix  G(t 1)  scaled 
by  hi.  Then,  only  the  mth  element  in  ( GH(fi)G{fi))~1GH{f1)H{l )  is  non-zero  (see  Eq.  (34)). 

If  the  Zth  channel  path  is  missampled,  i.e.,  77  f  fn,  'in  6  [1,  L\,  all  estimated  channel  taps 
{hn}n=\  are  generally  non-zero  even  if  noise  is  absent.  As  a  result,  energy  of  this  channel 
path  will  disperse  into  all  taps.  This  is  known  as  the  energy  leakage  phenomenon.  In  a 
multipath  channel  with  continuously  varying  path  delays,  the  equally-spaced  model  can 
not  simultaneously  sample  all  channel  paths.  Therefore,  energy  leakage  will  always  exist. 

The  tap  interval  Tp  is  set  as  the  inverse  of  the  signal  bandwidth,  which  is  known  as 
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the  time  resolution  of  the  system.  When  Tp  is  smaller,  GH(t\)G(f\ )  in  (33)  tends  to  be  ill- 
conditioned  and  the  problem  becomes  unsolvable. 

From  Eq.  (34),  the  nth  estimated  channel  tap  is  the  summation  of  contributions  from  all 
paths  of  the  physical  channel: 

L 

hn{fi)  =  Y.  hnjfi,  l )  +  77w(fi)  ,  n  =  1, . . . ,  L  (35) 

i=i 


where  hn(ri,  l)  contains  the  information  of  the  Ith  path  and  rjn ( T\ )  is  the  noise.  When  Tp  is 
chosen  as  the  inverse  of  the  signal  bandwidth,  hn(f\ ,  /)  can  be  further  expressed  by 


K(t i,l)  =hiex p 


x 


jn(N  -  1  )(rn  -  n) 
NTp 

sin(7r(fra  -  ti)/Tp) 
iYsin(7r(fn  -  ti)/(NTp)) 


(36) 


Eq.  (36)  is  actually  the  sampled  version  of  the  discrete  sine -function.  Therefore,  we  have 
the  following  results. 


Proposition  5  The  mth  channel  estimate  tap  contains  the  strongest  energy  from  the  Ith  channel 
path  if  m  =  argmin1<n<^  r;  —  rn\  ,  i.e.,  {rm,  hm}  is  the  closest  tap  to  {ti,  hi}.  Given  this,  energy 
leakage  on  the  other  taps  decreases  as  their  tap  indices  n  (n  f  m)  deviate  from  m. 


Proposition  6  Given  that  the  mth  estimated  channel  tap  contains  the  strongest  energy  from  {rg  hi}, 
hm(f\ .  l)  increases  as  \rt  —  fm\  decreases  in  [0,Tp/2).  Energy  on  the  other  taps  |^n(fi,/)|,  n  f  m 
approximately  decreases  when  \ n  —  fm\  decreases  in  [0,  Tp/ 2). 

Figs.  5(a)  and  5(b)  use  two  different  values  of  T\  to  channel  estimate.  We  can  see  that  by 
using  a  proper  f\,  energy  leakage  prior  to  the  first  path  {h\ ,  rx }  can  be  much  weaker. 


D.3  The  Proposed  ToA  Estimator 

For  the  ToA  estimation  purpose.  Fig.  5(b)  is  clearly  more  preferable  than  Fig.  5(a).  This  is 
because  when  the  channel  estimate  in  Fig.  5(a)  is  used  by  the  threshold-based  ToA  estimator, 
one  of  the  strong  energy  leakage  taps  can  be  mistakenly  picked  out  as  the  first  path.  This 
may  cause  a  severe  ToA  estimation  error.  Therefore,  we  need  to  suppress  the  energy  leakage 
before  estimating  the  ToA. 
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When  the  energy  leakage  is  sufficiently  suppressed,  a  sharp  jump  of  the  tap  amplitude 
will  emerge  near  the  leading  edge  of  the  channel  (see  Fig.  5(b)).  This  sharp  jump  of  am¬ 
plitude  can  be  detected  by  searching  the  value  of  77  to  maximize  the  following  energy  ratio 
between  two  adjacent  taps 

7n('Tl)  —  TT - -  .|9  ,n  e  [Li,L2\  (37) 

where  [L1;  L2]  represents  the  remaining  ToA  estimation  error  after  coarse  timing.  The  To  A 
estimate  is  then  obtained  by: 

A  =  T 1  +  (n  -  1  )TP  (38) 

with  (t\  ,  n)  =  argmax  7„(fi).  The  advantage  of  this  criterion  is  that  it  avoids  the 

0<ri <77,  L\<n<L2 

channel  dependent  threshold  required  by  the  threshold-based  ToA  estimators. 


D.3.1  Analysis  of  ToA  Estimation  Criterion 


Analysis  of  Eq.  (35)  for  an  arbitrary  multipath  channel  is  mathematically  intractable.  There¬ 
fore,  we  consider  a  simplified  two-path  channel  with  the  inter-path  interval  being  (p+0.5)Tp 
and  p  an  integer.  In  this  model,  the  first  path  carries  the  ToA  information  and  the  second 
path  models  the  interference  from  trailing  paths.  For  this  case,  the  strongest  interference 
arises  from  the  second  path  when  the  first  path  is  sampled  at  its  true  arrival  instant. 

Suppose  that  the  channel  paths  have  amplitudes  [hi,h2\  and  delays  [ti,t2].  Then  the 
energy  ratio  given  by  Eq.  (37)  can  be  expressed  as 


7n(Tl) 


IMri,  1)  +  hn(Ti,2)\2 

\hn-i(n,  1)  +  hn- i(ti,2)|2 


(39) 


where  hn(ri,  1)  and  hn(ri,  2)  are  contributions  from  the  two  channel  paths  [c.f.  (36)]. 

Using  t2  =  (p  +  0.5)TP  +  77  and  exp(j7r(p  +  0.5 )TP/ ( NTP ))  ~  1  when  p  <C  N,  we  have 

hn(fi,  2)  «  hn(n,  1)^7^ 

sin  (7r(fn  -  Ti)/(NTp))  sin  (7r(fn  -  t2)/Tp)  ^ 

sin  (7r(fn  -  n)/Tp)  sin  (7r(f„  -  t2)/(NTp)) 

Therefore,  energy  of  the  nth  tap  can  be  approximated  by  \hn(f\ )  |2  =  hn (t\  ,  1)| 2  +  \hn(fi,  2)|2 
and  the  energy  ratio  can  be  expressed  as 


7n(Tl) 


| M)l,  1)|2  +  \hn(Ti,  2)  | 2 

\hn-lin,  1)|2  +  \hn-l{fl,2)\2  ’ 


(41) 
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Since  Eq.  (41)  is  still  too  complicated  to  analyze,  we  need  to  further  simplify  it.  Suppose 
that  \hn(fi,  1)|  and  \hn(fi,2)\  reach  their  maximum  values  at  the  mth  and  m'th  taps  (m  < 
m'),  respectively,  such  that  \hm(fi,l)\  =  max(\hn(fi,  1)|)  and  |/w(t2,  1)|  =  max(|/v(A,  2)|), 
n,  n'  G  [1,  L\.  Numerical  analysis  has  confirmed  that  the  energy  ratio  will  only  be  maximized 
at  the  mth  and  m'th  taps  of  the  estimated  channel.  Based  on  this,  the  ToA  estimate  will  be 
either  fm  or  fm,  •  From  Proposition  5,  we  know  that  |ti  -  rm\  G  [0,  Tp/2)  and  |r2  -  rm>\  G 
[0,  Tp/2).  Certainly,  we  prefer  fm  than  fm' .  We  will  later  use  an  approach  to  avoid  the  energy 
ratio  being  maximized  at  the  m'th  tap.  For  now,  we  just  assume  that  the  ToA  estimate  is  fm 
and  analyze  yn(A)  at  the  mth  tap. 

Estimate  ToA  with  fm 

From  Proposition  6,  we  know  that  \hrn{f\ ,  1 ) |  is  a  decreasing  function  of  |ri  —  fm j.  Sim¬ 
ilarly,  \hm(f i,2)\  and  \hm-i  (fi  ,  2)  which  are  the  energy  captured  from  the  second  path  are 
increasing  functions  of  |r2  —  fm>\.  Due  to  the  relationship  r2  =  (p  +  0.5 )TP  +  T\,  |r2  —  fm>\ 
decreases  as  |ti  —  fm\  increases.  For  these  reasons,  both  \hm(fi,2)\  and  \hm_i(fi,2)\  are  in¬ 
creasing  functions  of  \hm(f\ ,  1)|.  From  (36),  the  maximum  value  of  \hm(f\,  1)|  is  hi  when  the 
first  path  is  exactly  sampled.  When  fm  =  (ti  —  0.5TP),  1)|  reaches  its  minimum  value 


\hm(ri,  1)| 


hi  sin(0.57r)  2  hi 

N  sin(0.57r/A)  n 


(42) 


In  order  to  relate  |/rm(ri,  1)|  with  |/im(ri,2)|  and  |/im-i(ri,  2)|,  we  use  the  following  ap¬ 
proximation 

2 


O  h- 

\hn(fi,  2) | 2  =  (  \hm(fi,  1)| - -  )  ce(n), 


IT 


(43) 


n  =  m,  m  -  1,  \hm{ru  1)|  G  [2Ai/7t,  hi 


where  2hi/n  is  the  minimum  value  of  \hm{fi,  1)|.  Coefficients  {ce(n)}((l=m-1  reflect  the  in¬ 
terference  strength  from  the  second  path.  In  particular,  the  stronger  the  second  path  or 
the  smaller  the  inter-path  interval,  the  larger  {ce(n)}™=m_1.  Numerical  analysis  shows  that 
ce(n)  is  almost  independent  with  \ti  —  fm\.  Therefore,  ce(n)  can  be  approximated  to  be  con¬ 
stant  to  |ri  —  fm |  and  accordingly  to  \hm(fi,  1)|.  In  addition,  numerical  result  shows  that 
ce(m  —  1)  ~  ce(m)  which  especially  holds  when  the  inter-path  interval  is  large. 

From  Proposition  6,  \hm-i{jii  1)1  decreases  as  \hm(f] ,  1)|  increases.  Based  on  this,  we  use 
the  approximation  \hm-i(fi,  1)|2  «  h\  —  hrn (t\  ,  1)|2.  Then,  the  energy  ratio  at  the  mth  tap  of 
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Figure  6:  (a)  ToA  estimation  by  the  maximizing  Eq.  (44).  (b)  ToA  estimation  by  maximizing  the 
energy  ratio  at  the  strongest  component  tap  of  the  first  path;  two  path  channel. 

the  estimated  channel  (41)  is  simplified  as 

\hm(n,  1)|2  +  (\hmin,  1)1  -  2hi/ir)2  ce{m) 

'in  (  \ 7)  « 

h\-\hrn(f1il)\2+(\hrn{fl,l)\-2h1/'K)  ce(m )  (44) 

1)|  G  [2/ii/ 7r ,  hi]. 

Eq.  (44)  contains  a  single  independent  variable  \hrn(f\ ,  1)|  with  the  coefficient  ce(m )  model¬ 
ing  the  interference  from  the  second  path.  For  each  ce(m),  the  ToA  estimate  can  be  calculated 
by  searching  \hm(fi,  1)|  to  maximize  the  energy  ratio  (44).  Fig.  6(a)  shows  that  the  ToA  es¬ 
timation  error  decreases  as  ce{m )  decreases,  which  is  a  quite  reasonable  result.  In  addition, 
when  ce(m )  is  very  small,  the  ToA  estimation  is  always  accurate.  When  ce(m )  is  large,  the 
ToA  estimation  error  will  rapidly  increase  because  the  second  path  is  dominant. 

Fig.  6(b)  shows  the  numerical  result  for  the  two-path  channel  where  the  ToA  is  esti¬ 
mated  by  maximizing  the  energy  ratio  of  (39)  at  the  mth  tap  that  contains  the  strongest 
contribution  from  the  first  path.  As  predicted  by  the  approximate  analysis  in  Eq.  (44),  the 
ToA  estimation  error  decreases  either  as  the  inter-path  interval  increases  or  the  second  path 
strength  decreases,  both  of  which  result  in  a  smaller  ce(m).  When  the  inter-path  interval  is 
large,  the  ToA  estimate  is  very  accurate  since  inference  from  the  second  path  is  very  weak. 
As  the  second  path  strength  increases,  the  ToA  estimation  error  will  rapidly  increase  since 
the  inter-path  inference  becomes  dominant. 

Avoid  the  Fake  ToA  Estimate  fm> 

In  the  preceding  analysis,  we  estimate  the  ToA  by  searching  f\  to  maximize  the  energy 
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Figure  7:  (a)  Avoid  picking  out  the  second  path;  M  =  4.  (b)  To  A  estimation  by  maximizing  the 
energy  ratio  among  the  reconstructed  taps  for  the  two  path  channel;  T2  =  n  +  3.5 Tp. 

ratio  7n(fi)  at  the  mth  tap  that  contains  the  strongest  contribution  from  the  first  path.  For 
this  approach  to  be  realistic,  we  need  to  avoid  the  energy  ratio  7n(fi)  being  maximized  at  the 
m'th  tap  that  contains  the  strongest  contribution  from  the  second  path.  We  use  the  following 
modified  energy  ratio 

7»(n, M)  = - hATJi - ,  „  6  [L,, L2\  .  (45) 

i—n—  1 

Different  from  (41),  the  denominator  in  (45)  also  includes  (M  —  1)  taps  prior  to  the  current 
tap.  Therefore,  if  M  is  sufficiently  large,  denominator  of  the  m'th  tap  will  not  only  include 
the  weak  energy  leakage  of  both  paths  but  also  the  strong  energy  leakage  from  the  first 
path  (see  Fig.  7(a)).  In  comparison,  denominator  of  the  mth  tap  only  contains  weak  energy 
leakage  from  both  paths.  As  a  result,  it  will  be  less  likely  for  Eq.  (45)  to  reach  its  maximum 
value  at  taps  that  contain  strong  energy  from  the  second  path. 

The  validity  of  this  method  has  been  confirmed  by  the  numerical  analysis  where  the 
To  A  is  estimated  by  searching  for  the  maximum  value  of  7n(fi,  M)  among  all  taps  of  the  es¬ 
timated  channel  (see  Fig.  7(b)).  In  Fig.  7(b),  an  optimal  value  of  M  exists  that  optimizes  the 
ToA  estimation  performance.  When  M  is  smaller  than  this  value,  ToA  estimation  accuracy 
improves  as  M  increases  since  it  becomes  less  likely  to  pick  out  the  trailing  paths  when  more 
leading  paths  are  included  in  the  denominator  of  7n(A,  M).  As  M  continuously  increases, 
ToA  estimation  performance  will  again  degrade.  This  is  because  as  more  weaker  energy 
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leakage  taps  prior  to  the  first  path  are  included,  the  estimator  becomes  less  sensitive  to  the 
ToA  estimation  error.  Simulations  show  that  these  results  also  hold  for  the  IEEE  802.15.4a 
UWB  channels. 

D.3.2  ToA  Estimation  for  Multipath  Channels 

Based  on  the  preceding  discussions,  we  have  the  following  ToA  estimation  algorithm  for 
general  multipath  channels. 

Proposed  Algorithm:  For  each  t1  e  [0 ,TP),  we  evaluate  the  energy  ratio  =  \hn(ji)\2 

/(if  1  \hi(fi)\2)  at  the  nth  tap,  n  e  [Li,  L2\.  [Li,L2\  represents  the  remaining  ToA  estimation 
error  after  the  coarse  timing.  Then  find  the  (77,  n)  pair  that  maximizes  the  energy  ratio  7n(fi,  M); 

(A,h)=  argmax  7„(fi  ,M).  (46) 

0<n<Tp,  Li<n<I/2, 

The  ToA  estimate  can  then  be  obtained  as: 

h  =  ri  +  {n  ~  1  )TP  .  (47) 

The  range  of  77  in  (46)  has  been  reduced  to  [0,  Tp)  as  compared  to  the  [0,  Th\  in  Eq.  (38)  where 
Th  is  the  maximum  channel  delay  spread.  This  is  because  for  the  equally-spaced  model,  the 
energy  ratio  satisfies  yn(A  +  kTp ,  M)  =  yn+k(^i,  M)  for  77  G  [0,  Tp)  with  k  being  an  integer 
(see  (36)).  Therefore,  it  is  sufficient  to  limit  77  in  [0,  Tp)  when  yn{f\ ,  M)  is  also  maximized 
with  respect  to  the  tap  index. 

D.3.3  FoFa  versus  Traditional  ToA  Estimators 

When  the  model-based  channel  estimator  is  used  for  ToA  estimation,  it  can  not  achieve  the 
best  accuracy  provided  by  the  channel  condition.  As  FoFa,  the  model-based  estimator  also 
uses  (33)  to  estimate  the  channel.  However,  it  does  not  consider  the  optimal  choice  of  77  in 
the  sense  of  ToA  estimation.  From  Fig.  5(a),  we  can  see  that  with  a  bad  choice  of  77,  delay 
of  any  of  the  first  5  taps  could  be  used  as  the  ToA  estimate,  due  to  their  sufficiently  large 
amplitudes.  In  such  cases,  the  resulting  ToA  estimation  error  will  be  comparable  to,  or  even 
much  greater  than,  the  tap  spacing.  Actually,  with  an  intention  to  solve  this  problem,  the 
FoFa  ToA  estimator  can  be  regarded  as  an  enhanced  sample-spaced  estimator. 

As  introduced  previously,  all  computations  in  the  FoFa  ToA  estimator  are  performed  in 
baseband  which  avoids  the  manipulation  of  analog  waveforms  as  in  correlator-based  esti- 
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mators.  In  addition,  FoFa  uses  the  channel  frequency  response  coefficients  for  ToA  estima¬ 
tion  which  are  already  available  for  equalization  purposes  in  FDE  systems  such  as  OFDM 
and  SC-FDE.  Therefore,  the  FoFa  ToA  estimator  can  achieve  a  much  lower  hardware  cost 
by  building  on  existing  FDE  structures  with  minimum  alteration. 

Traditionally,  ToA  estimation  can  also  be  carried  out  via  channel  estimation  which  aims 
at  maximizing  the  likelihood  or  minimizing  the  error  between  the  estimates  of  all  paths 
and  the  true  ones.  Since  the  contribution  from  the  first  path  is  small,  these  estimators  often 
encounter  poor  convergence  or  suffer  from  local  optimum  problems  when  being  used  as 
ToA  estimators.  The  FoFa  ToA  estimator  directly  addresses  the  delay  estimation  of  the  first 
channel  path.  Therefore,  FoFa  can  outperform  the  traditional  channel-estimation-based  ToA 
estimators  in  terms  of  accuracy.  This  can  be  verified  by  simulations  for  FoFa  and  the  well- 
known  SAGE  algorithm.  Furthermore,  the  FoFa  ToA  estimator  is  computationally  efficient 
since  we  focus  on  the  estimation  of  a  single  parameter. 

D.4  Summary 

In  this  research,  we  propose  a  high-precision  ToA  estimator  by  focusing  on  the  first  arrival 
(FoFa).  The  ToA  estimator  relies  on  a  novel  criterion;  that  is,  locating  the  first  channel  path 
by  minimizing  the  energy  leakage  prior  to  the  first  path.  By  directly  addressing  the  estima¬ 
tion  of  the  first  channel  path,  the  proposed  ToA  estimator  can  outperform  the  traditional 
SAGE  algorithm  in  terms  of  accuracy  and  computational  efficiency.  Furthermore,  the  pro¬ 
posed  ToA  estimator  operates  in  baseband  to  avoid  the  complicated  manipulation  of  the 
analog  or  oversampled  waveform  at  the  receiver,  and  can  be  implemented  with  minimum 
alteration  of  existing  FDEs. 

E  Ranging  with  Multi-Band  UWB  Signals:  Coherent  vs.  Non¬ 
coherent  Combination 

In  our  approach,  the  tap  spacing  Tp  dictates  the  resolution  of  the  candidate  set  At 

the  first  glance,  it  appears  that  the  timing  resolution  can  be  arbitrarily  improved  by  reducing 
Tp.  However,  when  Tp  is  smaller  than  the  inverse  of  the  bandwidth  occupied  by  the  channel 
frequency  response  coefficients,  the  matrix  GH{fi)G(fi)  in  (33)  tends  to  be  ill-conditioned 
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and  the  problem  becomes  unsolvable.  This  is  because,  even  though  the  Vandermonde  ma¬ 
trix  G{f\)  always  has  full  column  rank,  columns  of  G(f\ )  become  increasingly  correlated  as 
Tp  decreases.  In  the  extreme  case  when  Tp  — >  0,  the  rank  of  G(f i)  becomes  one.  For  this 
reason,  the  achievable  resolution  of  the  estimator  is  constrained  by  the  used  bandwidth. 
This  bandwidth  versus  resolution  problem  has  also  been  mentioned  in  the  literature  in  a 
different  context. 

For  multi-band  systems  such  as  the  MB-OFDM  UWB,  two  strategies  can  be  adopted  to 
combine  the  channel  information  from  subbands:  the  coherent  combining  and  the  noncoherent 
combining.  For  the  coherent  combining,  estimates  of  channel  frequency  response  coefficients 
for  all  subbands  are  jointly  used  to  estimate  the  time  domain  channel  with  Eq.  (33).  For 
the  noncoherent  combining,  the  time  domain  channel  is  first  estimated  for  each  subband. 
Then,  the  ToA  estimates  obtained  from  all  subbands  are  averaged  for  a  better  final  estimate. 
The  bandwidth  used  by  the  coherent  combining  is  larger  than  each  individual  subband.  In 
addition,  the  smallest  model  resolution  Tp  is  the  inverse  of  the  bandwidth  that  is  used  for  the 
time  domain  channel  estimation.  The  coherent  combining  can  provide  a  better  resolution 
than  the  noncoherent  combining. 

It  should  be  noted  that,  as  other  ToA  estimation  techniques  in  the  literature,  the  coherent 
combining  in  the  FoFa  ToA  estimator  also  requires  that  no  random  phase  rotation  exists  in 
subband  signals  after  carrier  demodulation.  If  there  is  random  phase  rotation,  one  would 
have  to  adopt  the  noncoherent  combining.  In  addition,  the  resolution  enhancement  by  the 
coherent  combining  is  obtained  at  the  price  of  increased  computational  complexity  by  com¬ 
puting  the  inverse  of  GH(fi)G(fi)  of  a  greater  dimension. 

We  simulate  the  performance  of  the  FoFa  ToA  estimator  based  on  the  UWB  MB-OFDM 
system  specified  in  the  ECMA-368  standard.  In  this  standard,  the  entire  UWB  spectrum  is 
divided  into  14  equally-sized  subbands.  For  each  subband  with  a  bandwidth  of  528  MFIz, 
the  multicarrier  modulation /demodulatioin  is  performed  with  a  128  point  IFFT/FFT.  A 
total  of  122  subcarriers  are  used  as  data,  guard  and  pilot  subcarriers.  Each  simulation  is 
carried  out  in  2000  randomly  generated  channel  realizations.  As  the  performance  of  the 
subspace-based  algorithms  are  usually  quite  sensitive  to  channel  order  mismatch,  we  will 
only  compare  the  performance  of  the  FoFa  ToA  estimator  with  the  SAGE  algorithm. 

We  first  compare  FoFa  with  the  SAGE  algorithm  in  the  IEEE  802.15.4a  LoS  office  channel 
(CM3)  (see  Fig.  9(a)).  Assume  that  the  remaining  timing  ambiguity  is  ±4.5  ns  after  the 
coarse  synchronization.  The  M  in  (45)  is  set  so  that  the  average  energy  is  calculated  in  the 
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Figure  8:  (a)  The  RMSE  performance  of  SAGE  and  the  FoFa  ToA  estimator  in  LoS  channels,  (b)  The 
RMSE  performance  of  SAGE  and  the  FoFa  ToA  estimator  in  NFoS  channels. 

range  of  2  ns.  For  a  fair  comparison,  the  channel  estimates  obtained  by  using  Eq.  (33)  with 
a  random  f\  are  fed  to  the  SAGE  ToA  estimator  as  the  initial  state  of  the  iterative  algorithm. 

Fig.  9(a)  shows  the  root-mean-square  error  (RMSE)  performance  of  FoFa  and  SAGE 
when  the  3  subbands  are  either  coherently  or  noncoherently  combined.  For  both  estimators, 
the  RMSE  is  much  smaller  than  the  receiver  sampling  interval  1.89  ns=l/528  MFIz  at  high 
SNR.  For  both  coherent  and  noncoherent  combining,  SAGE  outperforms  the  FoFa  estimator 
at  low  SNR.  This  is  because  the  SAGE  algorithm  estimates  the  delays  and  amplitudes  of  all 
channel  paths.  At  high  SNR,  FoFa  outperforms  SAGE  by  directly  estimating  the  ToA.  With 
our  setup,  the  Matlab  simulation  time  of  the  SAGE  estimator  is  about  10  times  that  of  FoFa1. 
Although  the  simulation  codes  are  not  optimized,  this  can  somehow  show  that  the  FoFa  ToA 
estimator  is  more  computationally  efficient.  In  Fig.  9(a),  we  also  compare  the  performance 
of  both  coherent  combining  and  noncoherent  combining.  As  predicted  in  the  preceding 
analysis,  the  coherent  combining  shows  a  better  resolution  (high  SNR  performance)  than 
the  noncoherent  combining  due  to  the  improved  resolution  of  {hi,  fi}f=1. 

The  FoFa  and  SAGE  ToA  estimators  are  then  evaluated  and  compared  in  the  IEEE  802.15.4a 
non-line-of-sight  (NLoS)  office  channel  (CM4)  (see  Fig.  9(b)).  Compared  to  the  LoS  channel 
environment,  both  estimators  show  performance  degradations.  Notice  that  this  is  also  a 
common  problem  for  existing  timing  and  synchronization  algorithms,  because  in  the  NLoS 
scenario,  the  first  channel  path  is  not  necessary  to  be  a  strong  one  and  the  synchronization 

1For  SAGE,  we  have  used  the  time  of  one  iteration  for  comparison. 
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Figure  9:  (a)  ToA  estimation  for  the  outdoor  LoS  and  NLoS  channels,  (b)  ToA  estimation  for  the 
industrial  LoS  and  NLoS  channels. 

is  more  likely  to  be  affected  by  noise  and  the  trailing  paths. 

We  further  test  our  algorithm  in  outdoor  and  industrial  UWB  channels.  Figs.  9(a)  and 
9(b)  show  that  performance  in  LoS  channels  is  better  than  the  corresponding  NLoS  chan¬ 
nels.  This  is  because  in  the  NLoS  channel,  the  first  path  may  not  be  strong  and  the  estimator 
is  more  affected  by  the  trailing  paths.  However,  there  is  an  exception  for  the  outdoor  envi¬ 
ronment  at  high  SNR.  Comparing  their  channel  realizations,  we  find  that  channel  paths  in 
CM6  are  weaker  but  sparser  than  CM5.  As  a  result,  the  CM6  channel  can  be  better  resolved 
than  CM5  and  the  localization  of  its  first  path  will  be  less  interfered  by  its  trailing  paths. 

The  proposed  ToA  estimator  is  also  compared  with  the  SAGE-based  ToA  estimator. 
Channel  estimate  results  obtained  by  (33)  with  a  random  A  are  fed  to  the  SAGE-based  ToA 
estimator  as  the  initial  state  of  the  iterative  algorithm.  Simulations  show  that  FoFa  out¬ 
performs  SAGE.  SAGE  is  a  practical  method  to  estimate  multiple  parameters  due  to  its  fast 
convergence  than  the  EM  algorithm.  However,  SAGE  has  the  local  optimum  problem,  espe¬ 
cially  when  the  number  of  unknown  parameters  is  large.  Therefore,  the  estimated  channel 
may  not  be  very  close  to  the  real  channel  at  the  first  path.  The  proposed  ToA  estimator 
directly  estimates  the  first  path  delay,  without  unnecessarily  caring  about  the  estimation 
errors  for  the  trailing  paths  and  therefore  enables  a  better  ToA  estimation  performance. 
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F  Ranging  with  Multi-Band  UWB  Signals:  Random  Phase 
Ratation 

Multi-band  orthogonal  frequency-division  multiplexing  (MB-OFDM)  UWB  systems  com¬ 
bine  OFDM  and  frequency  hopping  that  periodically  alters  the  carrier  frequency  among 
multiple  subbands.  Due  to  their  huge  bandwidth,  UWB  systems  can  realize  high  resolu¬ 
tion  time-of-arrival  (TOA)  estimation  for  time-based  wireless  ranging  and  localization.  For 
MB-OFDM,  this  capability  can  be  further  enhanced  if  channel  information  is  collected  from 
multiple  subbands  to  estimate  the  TOA. 

TOA  estimation  is  usually  carried  out  in  two  steps.  First,  the  time  domain  channel  is 
estimated.  Then,  the  first  path  is  detected  from  the  channel  estimate  and  its  delay  is  used  as 
the  TOA  estimate.  In  order  to  achieve  high  resolution  channel  estimation,  channel  informa¬ 
tion  from  multiple  subbands  should  be  coherently  combined  so  that  consecutive  subbands 
are  treated  as  a  single  larger  band.  Coherent  combining  requires  that  subband  signals  have 
the  same  phase  rotation  after  carrier  demodulation.  However,  this  cannot  be  guaranteed  by 
MB-OFDM  receivers  due  to  the  uncertain  initial  phase  states  of  modulator  and  demodulator 
oscillators.  Therefore,  it  is  necessary  to  develop  a  signal  processing  algorithm  to  calibrate 
random  phase  rotations  across  subbands. 

The  proposed  phase  rotation  calibration  algorithm  relies  on  the  energy  leakage  phe¬ 
nomenon  which  has  been  used  for  precise  TOA  estimation  as  discussed  in  the  preceding 
section.  Energy  leakage  is  essentially  due  to  the  limited  resolution  of  finite  signal  band¬ 
width  when  the  channel  is  estimated.  Analysis  in  this  work  indicates  that  in  the  presence 
of  random  phase  rotations,  resolution  of  MB-OFDM  signals  degrades  which  is  reflected  by 
increased  energy  leakage  in  the  channel  estimate.  Based  on  this,  we  propose  to  calibrate 
phase  rotations  of  subbands  by  suppressing  the  energy  leakage  effect.  Validity  of  the  pro¬ 
posed  technique  is  corroborated  in  IEEE  802.15.4a  UWB  channels. 

F.l  MB-OFDM  System  Model 

At  the  MB-OFDM  transmitter  (see  Fig.  10(a)),  information  bits  are  mapped  to  constellation 
points  and  multi-carrier  modulated  by  inverse  fast  Fourier  transform  (IFFT).  In  order  to 
avoid  the  inter-symbol  interference  (ISI),  guard  interval  (GI)  is  added  to  each  OFDM  sym¬ 
bol.  After  digital-to-analog  conversion  (DAC),  baseband  signals  are  carrier  modulated  and 
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transmitted  from  antenna.  At  the  receiver  (see  Fig.  10(b)),  received  signals  are  carrier  de¬ 
modulated  and  analog- to-digital  converted  (ADC).  After  coarse  timing,  GI  is  removed  and 
the  baseband  signals  are  multi-carrier  demodulated  by  FFT. 

MB-OFDM  systems  can  periodically  alter  the  carrier  frequency  according  to  the  time- 
frequency  code  known  to  both  transmitter  and  receiver  (see  Fig.  11).  Due  to  the  uncertain 
phases  of  transmitter  and  receiver  oscillators  (0t  and  <pr  in  Fig.  10),  baseband  received  sig¬ 
nals  contain  a  random  phase  rotation  (f>  =  [oL  —  <f>r)  which  varies  in  the  range  between  0 
and  2tt.  Although  the  random  phase  rotation  can  even  change  when  the  system  twice  en¬ 
ters  the  same  subband,  we  mainly  focus  on  the  inter-band  phase  calibration  across  different 
subbands.  After  that,  we  will  briefly  introduce  how  to  perform  the  intra-band  phase  cali¬ 
bration. 

Same  as  the  single-band  OFDM,  MB-OFDM  can  convert  the  frequency  selective  channel 
to  parallel  flat  fading  channels,  each  corresponding  to  a  subcarrier.  Suppose  that  the  system 
operates  over  B  subbands  and  each  subband  contains  K  subcarriers.  The  received  signal 
r^k  on  the  kth  (0  <  k  <  K  —  1)  subcarrier  of  the  6th  (0  <  b  <  B  —  1)  subband  can  be  expressed 
as: 

h>,fc  'A,fc-^f&,fcOXp  (j0ft)  -f-  ^b,k  •  (48) 

In  Eq.  (48),  s^k  is  the  signal  transmitted  on  the  /.  th  subcarrier  of  the  6th  subband  and  G./,-  is 
noise.  Hkj,  is  the  channel  Fourier  transform  coefficient: 

L—l 

Hb,k  =  ^  hiexply-jUb^ri)  (49) 

;=o 

where  {hi}^  and  {t/I/i,,1  are  multipath  amplitudes  and  delays  of  the  channel  impulse 
response  (CIR) 

L—l 

h(t)  =  ^2  -  n) ,  (50) 

1=0 

and  ub,k  is  the  frequency  of  the  (k,  6)th  subcarrier,  d/,  is  the  random  phase  rotation  in  base¬ 
band  signals  of  the  6th  subband. 

Based  on  (48),  channel  information  can  be  obtained  in  frequency  domain  by  the  follow¬ 
ing  operation: 

Hh,k  =  —  =  Hb,kexp(j(j)b)  +  rjb,k  (51) 

where  r]b,k  =  £b,k/sb,k  is  noise. 
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Figure  10:  MB-OFDM  system  model. 

F.2  Channel  Estimation  with  Coherent  Combining 

The  random  phase  rotation  problem  is  encountered  in  channel  estimation  with  coherent 
combining  which  is  the  first  step  of  high  resolution  TOA  estimation  for  MB-OFDM.  In  the 
absence  of  a  priori  channel  information,  the  maximum  likelihood  (ML)  criterion  results  in 
the  optimal  channel  estimator.  Given  that  in  Eq.  (48)  is  white  Gaussian,  the  ML  channel 
estimator  obtains  estimates  of  {/i;}[l01,  {ri  }!'=('  and  nuisance  parameters  {(j)b}b=o  when  the 
following  squared  error 

B- 1  K—l 

E  =  J2Y1  kfc-rVl2  (52) 

b= 0  k= 0 

is  minimized  between  received  signals  and  their  reconstructed  versions 

L— 1 

n,k  =  sb:k  T  hiexp(—jujbtkTi)exp  (j$b)  .  (53) 

;=o 

The  ML  estimator  is  too  computationally  intensive  to  be  feasible  for  UWB  channels  due  to 
the  huge  number  of  multipath  components. 

In  order  to  maintain  low  complexity. 


many  channel  estimation  related  problems 

■g 

are  investigated  based  on  the  sampled  chan-  J 
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nel.  This  suggests  that  the  time  domain  & 
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channel  can  be  reconstructed  by  fitting  fre¬ 
quency  domain  channel  information  (51)  to 
the  following  sequence  with  equally-spaced 
samples: 

L- 1 

h(t)  =  ^2  hn8(t  -  Tn)  •  (54) 

n= 0 

The  sequence  h(t)  contains  L  samples  with 

amplitudes  {hn}^Z g  and  delays  {fn}rnZ(\.  Given  the  first  sample  delay  f0  and  sample  spacing 
Tp/  delays  of  other  samples  can  be  determined  as  fn  =  f0  +  nTp,  1  <  n  <  L  —  1.  In  order  to 
reconstruct  the  whole  CIR,  LTP  should  be  larger  than  the  maximum  channel  delay  spread. 
In  this  work,  tap  spacing  Tp  is  chosen  as  the  inverse  of  the  entire  bandwidth  of  all  subbands 
which  is  known  as  the  system  resolution.  Compared  to  the  ML  criterion  (52)  and  (53),  this 
estimator  has  much  lower-complexity  because  it  has  only  one  free  parameter  f0.  In  the 
preceding  section,  f0  is  optimized  for  TOA  estimation  when  the  mistiming  induced  energy 
leakage  is  minimized. 

For  each  arbitrary  f0,  sample  amplitudes  of  sequence  (54)  can  be  obtained  by  fitting 
{hn}n=o  1°  the  frequency  domain  channel  information  (51)  as 


L—l 

5>„exp(  J^b^k^n)  -^b.kj 
n= 0 


(55) 


0  < k < K  —  h  0<6<5-l 


As  all  subbands  have  been  incorporated  into  channel  estimation,  Eq.  (55)  results  in  a  coher¬ 
ent  combining  of  MB  channel  information. 

The  relationship  (55)  can  be  expressed  in  matrix  form  as 


Gh  =  H 


(56) 


where  H  =  [H0j 0,  //o.i , . . . ,  Hb-i,i<-i]t  contains  all  frequency  domain  channel  information, 
h  =  [h0,  hi, ,  hi_ i]T  is  the  vector  of  sample  amplitudes  and  G  is  the  Fourier  transform 
matrix  determined  by  Eq.  (55).  A  simple  least  squares  (LS)  solution  of  (56)  can  be  formed  as 

h=  (GHG)~1  GUH  .  (57) 

It  should  be  noted  that  h  depends  on  the  first  sample  delay  r0  through  G. 
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Figure  12:  (a)  Amplitude  of  the  ith  path's  information  in  the  estimated  channel  in  the  presence  of 
random  phase  rotations,  (b)  Amplitude  of  the  ith  path's  information  in  the  estimated  channel  in  the 
absence  of  random  phase  rotations. 

As  mentioned  above,  sample  spacing  Tp  is  chosen  to  be  the  inverse  of  the  entire  band¬ 
width.  There  are  two  reasons  for  this.  First,  if  Tp  is  smaller  than  this  value,  GHG  in  Eq.  (57) 
tends  to  be  ill-conditioned  and  the  problem  becomes  unsolvable.  Secondly,  given  this  Tp, 
GhG  is  an  identity  matrix  and  the  computational  cost  is  reduced  by  avoiding  calculating 

{GnG)-\ 


F.3  Analysis  and  Calibration  of  Random  Phase  Rotation  Effect 

In  this  section,  we  first  analyze  the  impact  of  random  phase  rotations  to  channel  estima¬ 
tion.  Based  on  the  analysis  result,  we  will  propose  our  random  phase  rotation  calibration 
algorithm. 

F.3.1  With  Random  Phase  Rotations 

Because  Eqs.  (49)  and  (55)  are  linear  to  the  L  channel  multipath  components,  the  nth  sample 
hn  of  the  estimated  channel  h(t)  is  the  superposition  of  information  obtained  from  all  L 
paths 

L—l 

hn  =  ^2  hn(l)  +fjn,0<n<L  —  1.  (58) 

1=0 
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In  Eq.  (58),  rjn  is  the  noise  term  and  hn(l)  contains  channel  information  obtained  from  the 
/th  multipath  component  which  can  be  expressed  by 


B—l 


hn{l )  =  hi  y^exp(j0fc)exp 


6=0 


2nb 

]BTn 


n) 


x  hsn{l) 


(59) 


and 


KH) 


exp 


( J7T(/t  -  1  )(rw  -  T;)\ 

v  ^  / 


sin(7r(rn  -  ti)/(BTp )) 
iVsin(7r(rn  -  ti)/(NTp)) 


(60) 


In  Eq.  (59),  hn(l )  contains  channel  information  collected  from  all  B  subbands  distorted 
by  random  phase  terms  {exp{jOb)}l^('  with  {4>b}b=o  being  uniformly  distributed  in  the  rage 
between  0  and  2n.  Energy  captured  by  hn(l)  is  random  but  its  mean  value  can  be  obtained 
as  (cf.  (59)): 


e(|moT) 


h2  sin2(7 T(Tn-Ti)/(BTp)) 

1  KN sin2(7r(fn  -  ti)/(NTp )) 


(61) 


Note  that  the  energy  in  Eq.  (61)  is  determined  by  the  discrete  sine  function  sampled  at  time 
instants  (fn  —  p),  0  <  n  <  L  —  1. 

Fig.  12(a)  shows  an  example  of  Eq.  (61)  where  the  blue  line  represents  the  /th  channel 
path  with  delay  ti  =  0  and  red  lines  are  square  root  of  Eq.  (61)  with  their  envelope  (dashed 
line)  being  the  discrete  sine  function.  The  sine  function  has  a  main  lobe  width  of  2 BTP. 
From  the  definition  of  Tp,  BTP  is  the  inverse  of  bandwidth  of  a  single  subband.  This  implies 
that  no  matter  how  many  subands  are  used,  in  the  presence  of  random  phase  rotations, 
resolution  of  channel  estimation  is  only  facilitated  by  subband  bandwidth.  As  shown  in 
Fig.  12(a),  even  for  a  single  channel  path,  multiple  samples  in  the  estimated  channel  contain 
non-zero  energy.  This  is  known  as  the  energy  leakage  phenomenon.  It  should  be  noted  that 
due  to  the  random  phase  rotation  effect,  energy  leakage  cannot  be  avoided  by  horizontally 
shifting  the  sequence  (see  Fig.  12(a)). 


F.3.2  Without  Random  Phase  Rotations 

In  this  case,  we  assume  that  phase  rotations  of  all  subbands  have  been  perfectly  calibrated 
to  the  same  value,  for  example  0o-  Then  information  of  the  /th  path  in  the  estimated  channel 
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can  be  expressed  as  (cf.  (59)) 


7  m  i  { jn(N  -  l)(rn  -  Ti) 

hn{l )  =/i«exp(j0o)exp  ^ - — - 

sin(7r(fn  -  Ti)/Tp) 

X  N sm(n(fn  -  ti)/(NTp)) 


(62) 


Energy  captured  by  the  nth  sample  from  the  Zth  channel  path  is  also  a  sampled  discrete  sine 
function  but  with  a  smaller  main  lobe  width  of  2 Tp  (see  Fig.  12(b)): 


Tl  m|2  ,2  sm\n(Tn-Ti)/Tp) 

nU|  1  N2  sin2(7r(fn  —  ri)/(NTp))  ' 


(63) 


Same  as  (61),  Eq.  (63)  still  shows  energy  leakage  dispersed  from  the  Zth  path  over 
{ hn  ( Z ) } n=o  •  However,  when  the  Zth  path  is  exactly  sampled  by  the  sequence  if  a  proper 
r0  e  [0,  Tp)  is  chosen,  energy  leakage  in  Eq.  (63)  vanishes  and  all  energy  of  the  Zth  path  is 
captured  by  the  sample  that  has  delay  of  p.  This  is  not  a  property  possessed  by  (61)  when 
random  phase  rotations  exist. 

The  difference  between  (61)  and  (63)  indicates  that  energy  leakage  can  be  induced  by 
not  only  missampling  but  also  the  random  phase  rotation  effect.  In  the  presence  of  random 
phase  rotations,  main  lobe  of  (61)  is  B  times  wider  than  (63)  which  results  in  a  more  dis¬ 
persive  energy  leakage.  These  suggest  that  {( t>b}b=o  can  Be  calibrated  by  suppressing  the 
additional  energy  leakage  induced  by  random  phase  rotations. 


F.3.3  Phase  Rotation  Calibration  Algorithm 

Fig.  5(a)  shows  the  real-valued  UWB  channel  and  amplitude  of  the  channel  estimate  (58) 
using  the  first  two  subbands  of  the  ECMA-368  MB-OFDM  system.  For  multipath  channels 
with  continuously  varying  delays,  energy  leakage  always  exists  because  the  equally-spaced 
sequence  cannot  exactly  sample  all  channel  paths.  Among  all  samples  in  the  channel  es¬ 
timate,  those  prior  to  the  first  channel  path  can  well  reflect  the  extent  of  energy  leakage 
because  no  multipaths  fall  in  this  range.  We  have  seen  that,  TOA  is  estimated  by  choosing 
a  proper  f0  so  that  the  missampling  induced  energy  leakage  contained  by  these  samples  is 
minimized.  For  phase  rotation  calibration,  energy  leakage  is  induced  by  both  missampling 
and  random  phase  rotations.  Therefore,  energy  leakage  prior  to  the  first  path  can  be  min¬ 
imized  by  choosing  proper  values  for  both  f0  and  {0 & jo1-  Based  on  this,  random  phase 
rotations  can  be  calibrated  across  subbands  by  the  following  algorithm. 
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Proposed  Algorithm:  After  coarse  timing,  select  M  samples  before  the  first  channel  path 

that  only  contain  energy  leakage.  For  each  f0  e  [0,  Tp)  and  e  [0, 2n),  be  [1,  B  —  1],  calculate  the 
total  energy  leakage  of  these  M  samples 

M—l 

£{r0,  {fb}b= L1)  =  l^l2  •  (64) 

n=0 

Then  find  f0  and  {f>b}b=f  that  minimize  the  energy  leakage 

(To,  {}b}b=i)  =  argmin  £(f0,  {Mb= 11)  .  (65) 

0<ro<Tp,  0<rj>h<2n 

Correct  phase  rotations  by  multiplying  the  bth  subband's  received  signals  with  exp^—jfb)- 

F.3.4  Discussions 

We  have  introduced  the  inter-band  phase  calibration  of  different  subbands.  It  should  be 
noted  that  phase  rotation  also  changes  when  the  system  twice  enters  the  same  subband. 
Suppose  that  <j)b(i )  and  +  1)  are  phase  rotations  when  the  system  for  the  /  th  and  ( i  +  l)th 
times  enters  the  6th  subband  and  channel  does  not  change  in  this  period.  Calibration  of 
fb  (f)  and  f>b(i  +  1)  can  be  simply  realized  by  using  the  phase  of  the  following  correlation: 

I< 

nf,  %  + 1)  =  rb,k(i)rt,ki.i  + x)  •  (66) 

k= 1 

In  Eq.  (66),  the  same  set  of  training  data  are  transmitted  for  received  signals 

{rb,k(i)}k=o  an<4  {rb,k(i  +  l)}^1.  Combining  Eqs.  (66)  and  (65),  random  phase  rotations 
can  be  calibrated  for  the  general  case  with  arbitrary  time-frequency  hopping  patterns. 

In  Sections  F.2  and  F.3,  we  have  focused  on  the  coherent  combining  where  multiple  sub¬ 
bands  are  treated  as  a  single  larger  band  for  channel  and  TOA  estimation.  Analysis  indicates 
that  when  phase  rotations  of  subbands  are  calibrated,  coherent  combining  can  achieve  the 
resolution  facilitated  by  the  entire  bandwidth.  In  comparison,  noncoherent  combining  first 
estimates  the  CIR  and/ or  TOA  for  each  subband  and  then  averages  subband  estimates  for 
a  better  final  estimate.  Analysis  proves  that  noncoherent  combining  can  reduce  the  mistim¬ 
ing  probability  of  TOA  estimation  when  more  subbands  are  used.  Noncoherent  combining 
does  not  require  the  calibration  of  random  phase  rotations  and  its  computational  complex¬ 
ity  is  lower  than  coherent  combining.  However,  it  can  only  obtain  the  resolution  facilitated 
by  subband  bandwidth. 
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Figure  13:  (a)  Phase  calibration  performance  for  IEEE  802.15.4a  office  channels. (b)  Phase  calibration 
performance  for  IEEE  802.15.4a  industrial  channels. 

F.4  Simulations 

The  proposed  phase  rotation  calibration  algorithm  is  simulated  based  on  the  ECMA-368 
MB-OFDM  system  in  IEEE  802.15.4a  UWB  channels.  Channel  information  from  the  first  two 
subbands  are  coherently  combined  to  estimate  the  time  domain  channel.  M  —  4  samples 
are  used  to  calculate  the  energy  leakage  dispersed  prior  to  the  first  channel  path  (see  (64)). 

Figs.  13(a)  and  13(b)  show  performance  for  the  IEEE  802.15.4a  office  and  industrial  chan¬ 
nels  including  both  the  line-of-sight  (LOS)  and  non-line-of-sight  (NLOS)  scenarios.  Phase 
calibration  errors  normalized  by  27t  are  plotted  versus  the  signal- to-noise  ratio  (SNR).  For 
both  channel  types,  performance  in  the  LOS  scenario  is  better  than  the  corresponding  NLOS 
scenario.  This  is  because  that  in  LOS  channels,  the  first  channel  path  is  strong  and  energy 
leakage  (64)  is  dominated  by  the  first  path,  which  makes  LOS  channels  more  like  single 
path  channels.  Since  the  energy  leakage  based  algorithm  works  perfectly  for  single  path 
channels,  the  LOS  performance  should  be  better  than  NLOS. 

F.5  Summary 

In  this  work,  we  solved  the  random  phase  rotation  problem  which  is  encountered  in  the 
coherent  combining  of  MB  channel  information  for  high  resolution  TOA  estimation.  Our 
analysis  indicates  that  in  the  presence  of  random  phase  rotation,  resolution  of  MB  signals 
degrades  due  to  the  increased  energy  leakage  in  the  channel  estimate.  Based  on  this,  we 
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propose  to  calibrate  random  phase  rotations  of  subbands  by  minimizing  energy  leakage. 
Simulation  results  show  that  our  algorithm  works  well  in  IEEE  802.15.4a  UWB  channels. 


G  From  Ranging  to  Localization:  Cooperative  Localization 
via  SDP 

Recently,  GPS-complementary  location-aware  sensor  networks  have  attracted  increasing  in¬ 
terests  in  a  wide  range  of  applications  such  as  earthquake  detection,  weather  forecasting, 
current  flow  measurement.  Such  networks  typically  consist  of  two  types  of  sensors:  the 
anchors  or  base  stations  with  self-positioning  capability,  and  the  targets  with  their  positions 
unknown  to  the  system.  The  localization  problem  amounts  to  estimating  the  locations  of  the 
targets  based  on  the  known  positions  of  the  anchors  and  geometrical  relationships  among 
the  sensors. 

In  conventional  wireless  systems,  such  as  cellular  networks,  localization  is  achieved  in  a 
non-cooperative  manner;  that  is,  target  positions  are  estimated  only  by  the  anchors  within 
the  localization  range.  However,  sensor  networks  can  adopt  a  cooperative  mode  where 
inter-target  communications  are  also  permitted  such  that  targets  can  even  utilize  location 
information  of  the  anchors  outside  of  the  direct  communications  range  and  single-hop  con¬ 
nections  to  anchors  are  no  longer  a  must.  While  not  all  targets  can  estimate  their  locations 
non-cooperatively,  more  of  them  may  cooperatively  locate  themselves.  Cooperative  local¬ 
ization  not  only  reduces  the  density  of  the  expensive  anchors  but  also  provides  improved 
localization  accuracy.  Semidefinite  programming  (SDP)  technique  is  widely  adopted  to  con¬ 
vert  the  original  nonconvex  problems  into  convex  ones  in  recent  years.  Time-of-Arrival 
(ToA)  based  Standard  SDP  (SSDP)  was  proposed  in  the  literature.  However,  SSDP  is  not  in¬ 
herently  suitable  for  large-scale  networks  since  the  arithmetic  operation  complexity  of  SDP 
is  at  least  0(k3),  where  k  is  the  dimension  of  the  semidefinite  cones.  This  also  indicates  that 
the  semidefinite  cones  play  an  important  role  in  further  reducing  the  complexity  of  SDP.  In 
some  other  work.  Edge-based  SDP  (ESDP)  and  Node-based  SDP  (NSDP)  algorithms  were 
proposed  as  further  relaxations  of  SSDP  using  ToA  measurements.  This  method  of  relaxing 
a  set  of  low-dimensional  cones  instead  of  a  single  high-dimensional  cone  was  computation¬ 
ally  desirable.  Minimax  SDP  applicable  to  both  ToA  and  Received  Signal  Strength  (RSS) 
models  was  also  proposed. 
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Fundamental  limits  of  non-cooperative  and  cooperative  localization  with  ToA  measure¬ 
ments  were  exploited  recently.  In  this  research,  the  Cramer-Rao  Lower  Bound  (CRLB)  is 
tailored  for  our  SDP  algorithms  in  terms  of  the  system  models.  Besides  the  CRLB  for  ToA 
measurements,  we  also  derive  the  CRLB  for  RSS  measurements.  These  CRLBs  are  further 
used  to  show  that  cooperative  localization  has  the  potential  to  outperform  non-cooperative 
localization.  Our  minimax  Component-wise  SDP  (CSDP)  algorithm  is  applicable  to  both 
ToA  and  RSS  models  in  both  cooperative  and  non-cooperative  setups.  Here,  we  investigate 
the  effects  of  semidefinite  cones  in  terms  of  accuracy  and  efficiency  of  the  cooperative  lo¬ 
calization  algorithms  with  both  ToA  and  RSS  measurements.  Effects  of  SSDP,  ESDP,  NSDP 
and  CSDP  are  explored  and  analyzed  to  show  that  our  CSDP  is  the  best  among  them,  it 
is  efficient  while  retaining  the  crucial  theoretical  properties  of  the  SDP.  Simulations  will  be 
performed  to  verify  all  our  analyses. 

G.l  Cooperative  Localization  Algorithms 

In  a  sensor  network  localization  system  with  M  anchors  and  N  targets  in  a  D  dimensional 
real  Euclidean  space  1ZD  (D  =  2  or  3),  locations  of  the  anchors  constitute  a  known  vector 
set  Va  {o-i,  a2, . . . ,  aM},  and  locations  of  the  targets  form  an  unknown  vector  set  Vx  := 
{xi,  x2,  ■  ■  ■ ,  xN},  with  all  vectors  D-dimensional.  Furthermore,  we  are  given  the  Euclidean 
distances  dmn  between  am  and  x„  for  some  (m,  n),  and  d,3  between  xt  and  x:j  for  some 
Specifically,  let  £a  :=  {V  (m,n)  :  am  G  Va,xn  G  Vx,  dmn  is  specified},  £x  :=  {V  (i,j)  : 
Xi,Xj  G  14,  dtj  is  specified  and  i  <  j}.  The  localization  problem  in  1ZD  for  the  undirected 
graph  Q  :=  (V,  £),  where  V  :=  Va  U  Vx,  £  :=  £a  U  £xr  is  to  determine  the  coordinates  of  the 
unknown  target  positions  Vx  from  the  known  anchor  positions  Va  and  the  partial  distance 
measurements  £. 

Notation:  |  •  |  denotes  the  absolute  value,  ||  •  ||  denotes  the  I2  vector  norm,  tr{-}  is  the 
trace  of  a  square  matrix,  (-)T  is  the  matrix  transpose  operator,  and  diag(-)  represents  a  di¬ 
agonal  matrix.  Boldface  lowercase  letters  denote  vectors,  boldface  uppercase  letters  denote 
matrices.  Specifically,  0  denotes  an  all  zero  entry  vector,  e,  denotes  an  N  dimensional  vector 
with  all  of  its  entries  zero  except  a  unit  entry  at  the  i-th  position,  Idxd  denotes  a  D  x  D- 
dimensional  identity  matrix.  Vector  and  matrix  dimensions  will  be  clear  in  the  context,  and 
will  be  specified  whenever  necessary.  Z  A  0  means  that  Z  is  positive  semidefinite  (PSD), 
where  Z  is  a  square  matrix.  Zult denotes  the  principal  submatrix  of  Z  from  its  rows 
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and  columns  indexed  by  ii, . . . ,  ik.  \A\  denotes  the  number  of  elements  in  the  set  A.  A\B 
means  the  set  difference  {x\xi  e  dflSj  ^  B}. 

Distances  are  measured  so  long  as  the  node  pairs  are  within  the  communications  range 
of  each  other.  Suppose  there  are  K  available  distances,  out  of  which  K  are  between  one 
target  and  one  anchor,  and  the  remaining  are  between  two  targets. 


G.1.1  ToA 


Without  loss  of  generality,  we  can  denote  the  target-anchor  distances  as 

dmn  —  ||®m  Xn\\  T  ITlTOn  •  V (  ffl  .  Tl)  £  £a, 


(67) 


and  represent  the  target-target  distances  as 

<^j  =  || Xi  -  Xj\\2  +  e £  £x  and  i  <  j,  (68) 

where  noises  mmn  and  il,?  are  i.i.d.  Gaussian  distributed  with  zero  mean  and  variance  a2. 

Denote  the  residue  vector  as  e  =  [eq,  e2, . . . ,  £r,  D?+ i,  £r+ 2,  ■■■Ak\  =  [£mn,  ■  ■  ■  Aij,  ■  ■  ■]  = 
[d2mn-\\am-xn\\2 1 . . . ,  dij  —  Wxi—XjW2, . . .],  and  XDxN  =  [xi,x2,-..,  a;  at],  then  the  conditional 
probability  density  function  of  the  measured  distances  is 


^V(m,n)G^a  and  i<j  £ ij 

—  _ _ _ p  o  cr^ 

”  (2na2)K/2^  5 

where  Yl\/(m  n)&ea  £mn  represents  the  traditional  non-cooperative  information,  and 
Evfijjef,  and  i<j  £ij  represents  the  cooperative  information. 

The  Maximum  Likelihood  (ML)  estimate  of  X  is 


(69) 


X  =  argrmii  ^  s2mn  +  ^  e\  =  arg  nun  ||e||2.  (70) 

V(m,n)eA/’a  V(iJ)eA4  and  i<j 

Note  that  X  is  optimal  in  the  ML  sense. 

The  optimization  problem  of  (70)  is  highly  nonlinear  and  nonconvex  in  X,  thus  difficult 
to  solve.  In  this  work,  we  use  minimax  relaxation  and  semidefinite  relaxation  to  convert  the 
original  nonconvex  problems  into  convex  ones.  Denote  Y  =  X1  X,  and  approximate  this 
nonlinear  and  nonconvex  equation  via  semidefinite  relaxation  Y  'p  X1  X.  The  result  is  still 
nonlinear  but  already  convex.  The  minimax  SDP  algorithm  for  ToA  measurements  2  is 


min  t 

t,Z 


(71) 


2Suppose  all  targets  are  localizable.  Localizability  itself  is  a  challenging  research  area  which  is  beyond  the 
scope  of  this  work. 
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V  (m,  n)  6  4iV  (i,  j)  G  £x  and  i  <  j, 
Z(1  :D)  =  I DxDi  Z(1  :D,i)  h  0,  V  Xi  G  Vx\X , 
z(i :D,i,j)  h  0,  V  (i,j)  G  £x  and  i  <  j. 


G.1.2  RSS 

In  the  ideal  case,  each  target  node  emits  a  signal  with  power  P,  and  anchor  nodes  within 
the  localization  range  can  receive  a  signal  with  energy  strength 

Smn  =  -PH^m  ®n||  V(?7T,  Tl)  G  £a ,  (72) 

and  target  nodes  within  the  localization  range  can  receive  a  signal  with  energy  strength 

Sij  =  P\\xi  -  G  £x  and  i  <  j ,  (73) 

where  /3  is  the  path  loss  coefficient. 

Under  the  lognormal  fading  channel,  the  RSS  in  dB  can  be  modeled  as 

10  log  smn  =  10  log  P  -  10/3 log(||om  -  xn\\)  +  emn,  (74) 

V(m,  n)  G  £a  for  target-anchor  connections  and 

10 log =  10 log P  -  10/3 log(||®j  -  xj ||)  +  fi,j,  (75) 

V(i,  j)  G  £x  and  i  <  j  for  target- target  connections,  where  mmn  and  fil3  are  assumed  to  be 
i.i.d.  Gaussian  distributed  noises  with  zero  mean  and  variance  a2. 

Denote  the  residue  vector  as  e  =  [ei,  e2, . . . ,  er<,  eR+1,  eK+2,  ...,eK\  =  [emn, . . .]  = 
[In  smn  —  in(||Q  Px  p), ....  In  sVJ  —  1 1 i (  x  l’x — ^ ), . . .],  the  conditional  probability  density  func¬ 
tion  of  the  RSS  in  dB  is 


p(lnsmn, . . In  Sij,  ■  ■  •  l-X:) 

_  Sy(min)gga  emn+^V(iJ)  ggj;  cij 

=  1  r  2(^W}2 

(27r(^g^)2)K/2  ’ 

where  n)e£a  emn  represents  the  traditional  non-cooperative  information,  and 
Ev(i,j)  6£x  and  i<j  4  represents  the  cooperative  information. 


(76) 
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The  ML  estimate  of  X  can  be  obtained  as 


A  =  arg min  ^  e2mn  +  ^  =  arg  nun  ||e|j2.  (77) 

V(m.,n)eA/’a  V(ij')eA4 

Note  that  X  is  optimal  in  the  ML  sense. 

Again,  the  optimization  problem  of  (77)  is  highly  nonlinear  and  nonconvex  in  X,  thus 
difficult  to  solve.  Similar  to  the  ToA  case,  we  adopt  an  minimax  SDP  approach  to  convert 
the  original  nonconvex  problems  into  convex  ones.  The  minimax  SDP  algorithm  for  RSS 
measurements  is 


min  t 

t,Z 


V  (m,  n)  G  £a ,  V  ( i ,  j)  G  £x  and  i  <  j, 
Z(1  :D)  =  IdxD,  Z(l:D,i)  h  0,  V  Xi  G  Va;\A, 
Z(1  :D,i,j)  h  o,  V  (i,j)  G  £x  and  i  <  j. 


(78) 


In  (71)  and  (78)  A  =  {xj  \\/i,  3(i,  j)  G  £x,  i  <  j  U  (n,  i)  G  n  <  i}  is  the  set  formed  by  all  the 
targets  that  have  other  target(s)  as  neighbor(s).  The  first  conic  constraint  for  those  x%  G  X 
is  redundant,  since  it  is  implied  by  the  second  conic  constraint.  Once  the  solutions  of  (71) 
and  (78)  are  found,  the  upper  right  D  x  N  corners  of  Z  are  taken  out  as  the  semidefinite 
approximated  location  estimate  X. 


G.2  Cramer-Rao  Lower  Bound 

Here,  we  derive  the  Cramer-Rao  Lower  Bound  (CRLB)  for  both  cooperative  and  non- 
cooperative  scenarios  with  both  ToA  and  RSS  measurements  according  to  our  system  mod¬ 
els.  These  CRLBs  are  then  used  to  make  comparisons  between  cooperative  localization  and 
non-cooperative  localization. 
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G.2.1  ToA 


Let  %toA  =  [d2mni  •  •  • ,  fiy- , . .  .]7  be  a  vector  that  contains  all  measurements.  From  (69),  we  see 
that  %toA  is  Gaussian  distributed  with  mean  nToA  and  covariance  matrix  CtoA 


l^ToA  —  [II  am  xn\\  j  •  •  •  >  ||  xi  xj  1 1'  j  •  ■  •]  j 


CtoA  =  diag(cr2, . . . ,  a2, . . .). 


The  Fisher  information  matrix  (FIM)  has  elements 


tToA  = 


d^ToA  r'  1  dfiToA 
dXi  '-'ToA  dXj 


+  \^{CtIa 


dCToA  /">— 1 
a®,  '-'ToA 


dC  ToA 

dXj 


}, 


i,j  =  1,2,...,  AT. 


(80) 


Considering  given  assumptions  in  our  ToA  model,  the  FIM  is 

JToA  =  4cj_2IToj4.  (81) 

Diagonal  terms  of  ITOj4  are 

IT  ToA  _  V-  ,/2  ,  J2 

Tin  /-ora,V(m,n)G£a  L'”mn  '  Z^/itW(i,n)££x  and  i<n  “in 

"F  and  n<j  ^  1)2,...,  A,  (82) 

and  off-diagonal  terms  are 

1^°'4  =  =  -dfj,  V(f,  j)  G  and  <  <  j.  (83) 

Diagonal  terms  are  contributed  by  all  connections,  including  both  the  cooperative  and 
the  non-cooperative  connections,  while  off-diagonal  terms  display  the  effects  of  the  coop¬ 
erative  links.  Thus  elements  of  IToA  for  the  non-cooperative  case  with  ToA  measurements 
only  contains  non-zero  diagonal  terms 

E  ^n,n=l,2,...,iv,  (84) 

m,V(m,n)g£a 

with  all  off-diagonal  terms  zero. 

Suppose  every  targets  is  connected  to  at  least  one  anchor,  then  lToA  is  a  strictly  diagonally 
dominant  matrix.  By  Levy-Desplanques  theorem  ,  IToA  is  non-singular.  3  Thus,  the  CRLB 
for  X  with  ToA  measurements  can  be  obtained  from  the  diagonal  elements  of  the  inverse 
matrix  of  IToA. 

3If  a  target  has  no  target-anchor  connections,  then  without  cooperative  information,  it  will  be  an  isolated 
non-localizable  node.  Here  we  do  not  consider  this  type  of  cooperative  vs.  non-cooperative  comparisons. 
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G.2.2  RSS 


Denote  %rss  =  [lnsmn, ....  In  slJ, . .  .]T.  From  (76),  tlRSs  is  Gaussian  distributed  with  mean 
Prss  =  tIn(  \\am-xn\\p)'  ■■■■  ln(pA^)>  •  •  ']T  and  covariance  matrix  CRSS  =  diag((^)2, 
. . . ,  (^yy^)2,  •  •  •)•  Similarly,  the  FIM  is 


J-RSS  _  /  IQfi  \2  -2yRSS 
1  111  10  ’ 


(85) 


Where 


ttRSS  _  y-  1 _ i_  x  '  J_ 

nn  £-sm,V(rn,n)e.£a  _r  Z-/i,V(i,ri)e£x  andi<ra  ^ 

+  Y. <j,V(n,j)e£x  and  n<j  ~cp~.  '  n  ~  1;  2,  ....  A'.  (86) 

nj 

and 

=  ijE  =  -  4“,  V(i,  j)  G  and  i  <  j.  (87) 

% 

Similarly,  elements  of  IRSS  for  the  non-cooperative  case  with  RSS  measurements  only 
involves  non-zero  diagonal  terms 

Inn S  =  E  4-,n  =  l,2,...,N,  (88) 

m,V(m,n)s£a  mn 

with  all  off-diagonal  terms  zero. 

We  can  then  get  the  CRLB  for  X  with  RSS  measurements  from  the  inverse  matrix  of  the 
FIM  IRSS. 


G.2.3  Cooperative  vs.  Non-cooperative  Localization 

Based  on  the  CRLB  derived  in  Section  G.2.1  and  G.2.2,  we  show  that  performance  bounds 
of  cooperative  algorithms  will  never  be  above  the  non-cooperative  counterparts. 


Theorem  1  The  Cramer-Rao  Lower  Bound  (CRLB)  of  the  cooperative  localization  is  lower  than  or 
equal  to  that  of  the  non-cooperative  localization. 


CRLBTf0%  <  CRLBTCcoop 


CRLB 


RSS 


Coop  - 


<  CRLB 


RSS 

Noncoop  • 


(89) 


Proof:  Denote  eigenvalues  of  IF  oA  as  Ai,  A2, . . . ,  A at.  By  Gersgorin  disc  theorem , 


|An-IL 


To  A  | 
nn  I 


sE 

i^n 


tToA.  I 

ni  I 


,n  =  1,2, 


,N. 


(90) 
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For  those  xn's  that  have  no  cooperative  connections 


An  =  Y,  d 

m,V(ra,n)E£< 


2 

mn 5 


and  for  those  xn's  that  have  cooperative  connections 


^ n  ^  ^  ^  d> 

m,V(m,n)6£, 


2 

mn • 


Hence,  eigenvalues  of  the  inverse  matrix  of  VoA  are 

1 


A-1  = 


E 


rl  2 

m,V(m,n)££a  umn 


for  those  xn's  that  have  no  cooperative  connections  and 

A”1  < 


E 


d2 

m,V(m,n)££a  umn 


for  those  xn's  that  have  cooperative  connections. 
From  (93)  and  (94),  we  conclude  that 


CRLBg£p  <  CRLB^toop. 


Similarly,  for  RSS  measurements 


(91) 


(92) 


(93) 


(94) 


(95) 


CRLBgSp  <  CRLB^s0snc00p. 


(96) 


For  both  ToA  and  RSS  measurements,  we  get  the  conclusion  that  cooperative  algorithms 
have  the  potential  to  outperform  non-cooperative  ones. 


G.3  Effects  of  Semidefinite  Cones 

Since  objective  functions  and  all  inequality  constraints  of  various  minimax  SDP  algorithms 
remain  the  same,  to  compare  minimax  SDP  algorithms,  we  only  need  to  compare  all  the  ex¬ 
isting  and  our  proposed  semidefinite  cones.  We  then  investigate  into  effects  of  semidefinite 
cones  on  efficiency  and  accuracy  of  the  cooperative  localization  algorithms. 

G.3.1  Types  of  Semidefinite  Cones 

Standard  Semidefinite  Cones 
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Standard  semidefinite  approximation  is  Y  =  X1  X  =$■  Y  F  X1  X,  i.e.,  the  standard 
semidefinite  cone  is 

IdxD  XdxN 
X JfxD  ^iVxJV 

Further  Relaxations  to  Standard  Semidefinite  Cones 

•  Edge-based  Semidefinite  Cones 

z(i:D,i,j)  h  0,  V  (i,j)  g  Sx.  (98) 

•  Node-based  Semidefinite  Cones 

z(i:D,iMXi)  b  o,  V  Xi  g  Vx.  (99) 

where  A fXi  =  {j  ■  (■ i,j )  e  Sx,i  <  j  or  (j,i)  e  j  <  ?'}  is  the  set  formed  by  all  the 
targets  connected  to  xt. 

•  Component-wise  Semidefinite  Cones  (our  proposed) 

Z (1  :D,i)  —  O’  V  Xi  G  Vx\X . 

Z(l:D,i,j)  h  0,  V  (i,j)  G  £x-  (100) 


G.3.2  Performance  comparisons 

In  the  literature,  it  has  been  shown  that  FSSDP  C  Fnsdp  C  Fesdp,  where  F  represents 
the  solution  set  of  the  corresponding  SDP  algorithm.  For  our  proposed  CSDP,  we  have  the 
following  result. 

Proposition  7  The  solution  sets  ofNSDP,  CSDP  and  ESDP  follow  the  following  relationship: 

FNSDP  c  fCSDP  c  pESDP _  (101) 

Proof  Denote  solutions  to  (99),  (100)  and  (98)  as  Z*NSDP,  Z*CSDP  and  Z*ESDP,  respectively. 
NSDP  contains  redundancy  because  it  assumes  that  all  pairwise  connections  of  Afxt  exist, 
which  is  not  necessarily  true.  After  removing  these  unspecified  variables  in  Z*NSDP,  one 
would  obtain  Z*CSDP.  Similarly,  after  removing  the  non-cooperative  constraints  in  Z*CSDPr 
one  would  obtain  Z*ESDP. 

To  further  explore  the  performances  of  these  SDP  algorithms,  let  us  first  introduce  some 
definitions. 
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Definition  1  An  undirected  graph  is  a  chordal  graph  if  every  cycle  of  length  greater  than  three  has 
a  chord.  An  equivalent  definition  is  that  any  chordless  cycles  have  at  most  three  nodes. 

Definition  2  A  square  matrix,  possibly  containing  some  unspecified  entries,  is  called  partial  sym¬ 
metric  if  whenever  the  (■ i,j)-th  entry  of  the  matrix  is  specified,  then  so  is  the  ( j,i)-th  entry,  and 
the  two  are  equal.  A  partial  semidefinite  matrix  is  a  partial  symmetric  matrix  for  which  every  fully 
specified  principal  submatrix  is  positive  semidefinite. 

Given  these,  we  can  now  introduce  the  following  Lemma. 

Lemma  2  Every  partial  positive  semidefinite  matrix  with  undirected  graph  Q  has  positive  semidefi¬ 
nite  completion  if  and  only  if  Q  is  chordal. 

Based  on  the  above  proposition,  definitions  and  lemma,  we  have  the  following  theorem. 
Theorem  2  If  the  undirected  graph  Qx  :=  (Vx,  Sx)  induced  by  all  the  targets  is  chordal,  then 


pSSDP  _  pCSDP 


(102) 


Proof:  By  Proposition  7,  we  only  need  to  prove  that  any  solution  to  (100)  can  be  completed 
to  a  solution  of  (97).  Suppose 


Z  = 


(103) 


is  a  solution  of  (100).  By  Definition  2,  Z  is  a  partial  semidefinite  matrix. 

By  the  prerequisite,  the  undirected  graph  induced  by  Y  is  chordal.  We  now  prove  that 
the  undirected  graph  induced  by  Z  is  also  chordal.  Notice  that  the  graph  induced  by  Z 
has  N  +  D  nodes,  with  every  specified  off-diagonal  entry  denoting  an  edge.  Every  one  of 
the  first  D  nodes  in  the  graph  induced  by  Z  is  connected  to  every  target  node  in  Vx.  If 
a  cycle  in  the  graph  induced  by  Z  contains  some  or  all  of  the  first  D  nodes,  then  it  must 
have  a  chord  since  each  of  the  first  D  nodes  is  connected  to  every  target;  if  the  cycle  does 
not  contain  any  of  the  first  D  nodes,  then  it  must  still  contain  a  chord  because  the  graph 
induced  by  Y  is  chordal.  Thus,  the  graph  induced  by  Z  is  chordal.  By  Lemma  2,  Z  must 
have  a  positive  semidefinite  completion  Z.  Considering  that  (97)  and  (100)  share  the  same 
constraints  specified  by  the  existing  connections,  Z  must  be  a  solution  to  (97). 

For  chordal  graphs  Qx,  we  do  not  lose  any  relevant  information  in  CSDP,  thus  SSDP  and 
CSDP  are  equivalent.  Notice  that  we  only  need  Qx  to  be  chordal,  not  requiring  that  the 
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whole  localization  graph  Q  to  be  chordal.  Furthermore,  by  Proposition  7  and  Theorem  2, 
one  can  easily  deduce  that 

FSSDP  =  pNSDP  (104) 

for  a  chordal  graph  Qx,  which  coincides  with  the  results  in  the  existing  literature. 

These  result  essentially  confirm  that,  containing  reduced  constraints,  CSDP  leads  to 
identical  solutions  as  SSDP  and  NSDP  In  other  words,  both  SSDP  and  NSDP  contain  redun¬ 
dant  variables,  and  these  redundant  elements  increase  the  complexity,  thus  have  negative 
impacts  on  the  localization  efficiency.  ESDP  is  theoretically  deficient,  because  the  removed 
D  +  1-dimensional  conic  constraints  in  (98)  in  the  proof  of  Proposition  7  are  necessary  infor¬ 
mation.  ESDP  only  retains  cooperative  information  while  loses  non-cooperative  informa¬ 
tion. 

G.3.3  Complexity  comparisons 

In  Table  1,  we  list  the  computational  complexity  of  various  SDP  localization  algorithms. 
Note  that,  ^=1(D  +  l  +  \J\fXi I)2,  \£x\(D+2)2,  and  (N  -\X\)(D  +  l)2  +  \£X\(D  +  2)2  are  typically 
much  smaller  than  (N  +  D)2.  Hence,  with  a  much  smaller  number  of  variables,  NSDP,  ESDP 
and  CSDP  have  the  potential  to  reduce  the  computational  complexity  compared  to  SSDP, 
especially  for  large  N. 


Table  V.  SDP  localization  algorithm  complexity  analysis 


Number  of 

Conic  Constraints 

Dimension 

of  Cones 

Number  of 

Variables 

Number  of 

Equality  Constraints 

SSDP 

1 

N  +  D 

(N  +  D)2 

\sa\  +  \sx\ 

NSDP 

N 

D  +  1  +  |Aav 

zL(d  +  i  +  \MXi\y 

\Sa\  +  \Sx\ 

ESDP 

\sx\ 

D  +  2 

\Sx\(D  +  2)2 

\Sa\  +  \Sx\ 

CSDP 

N-  \X\  +  \Sx\ 

D  +  1,  D  +  2 

(N  —  \X\)(D  +  l)2  +  \£X\(D  +  2)2 

\Sa\  +  \Sx\ 

Generally  speaking,  ESDP,  NSDP  and  CSDP  are  all  further  relaxations  of  SSDP.  Even 
though  SSDP  is  more  accurate  and  may  possess  a  better  efficacy  for  small-scale  networks,  it 
is  not  suitable  for  large  scale  networks  since  it  involves  a  single  high-dimensional  semidef- 
inite  cone.  ESDP  loses  some  useful  information,  thus  it  is  not  robust  enough  and  blows 
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up  under  severe  noisy  cases.  Even  though  NSDP  involves  multiple  low-dimensional  conic 
constraints,  it  contains  redundancy.  CSDP  naturally  results  from  the  component-wise  for¬ 
mulation  of  the  original  problem,  keeps  all  the  specified  and  necessary  variables,  and  only 
those  variables.  Thus  CSDP  can  achieve  comparable  accuracy  with  reasonable  complexity. 
Moreover,  given  the  complexity  analysis,  CSDP  is  suitable  for  large-scale  wireless  sensor 
networks. 

G.4  Simulations 

In  this  section,  we  use  the  CVX  package  as  our  SDP  solvers.  Localization  accuracy  is  eval¬ 
uated  by  the  mean  square  error  (MSE  =  A  J2n=i  II ~  xn\\2)  of  the  location  estimates,  and 
localization  efficiency  is  indicated  by  the  average  time  used  for  locating  all  the  targets  in 
the  sensor  networks.  System  parameters  for  RSS  measurements  are  P  =  1000,  (5  =  3.  All 
simulation  results  are  averaged  over  2000  Monte  Carlo  tests. 

On  a  2-D  map,  anchors  are  located  at  Va  =  {[0, 6]T,  [0, 0]r,  [2, 4]T,  [6, 0]r,  [7,  7]T,  [8, 0]T}, 
targets  are  placed  at  14  =  {[0, 4]T,  [3, 8]r,  [2, 2]T,  [6, 5]T,  [7, 3]T}.  Communications  are  allowed 
between  node  pairs  at  £a  =  {(1,1),  (2,1),  (1,2),  (3,1),  (2,3),  (3,2),  (3,  3),  (5, 2),  (3,4),  (4,3), 
(5, 4),  (4,  5),  (6,  5)}  and  £x  =  {(4,  5)}.  Cases,  including  inside  the  convex  hull,  on  the  edge  of 
the  convex  hull,  out  of  the  convex  hull,  as  well  as  the  cooperation  between  two  targets,  are 
all  considered  in  this  localization  system.  We  use  the  cooperative  algorithms  to  make  com¬ 
parisons,  and  give  the  CRLBs  for  both  cooperative  and  non-cooperative  localization  show¬ 
ing  that  the  cooperative  localization  has  the  potential  to  outperform  the  non-cooperative 
counterpart. 

See  Figure  14(a)  for  algorithm  accuracy  performances  and  Figure  14(b)  for  algorithm 
efficiency  performances  with  ToA  measurements.  See  Figure  15(a)  for  algorithm  accuracy 
performances  and  Figure  15(b)  for  algorithm  efficiency  performances  with  RSS  measure¬ 
ments.  CRLBs  are  also  given  as  benchmarks,  showing  that  accuracy  performances  of  our 
SDP  algorithms  are  satisfactory.  Our  simulation  results  show  that  CSDP  achieves  compara¬ 
ble  accuracy  with  reasonable  complexity.  Moreover,  by  adopting  cooperative  localization, 
we  successfully  locate  all  the  targets,  including  x4  and  .x\5,  which  are  not  localizable  without 
cooperative  information. 
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MSE(m  )  MSE(m  ) 


(a)  (b) 

Figure  14:  (a)  ToA  localization  accuracy,  (b)  ToA  localization  efficiency. 


(a)  (b) 

Figure  15:  (a)  RSS  localization  accuracy,  (b)  RSS  localization  efficiency. 
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G.5  Summary 

In  this  research,  Time-of-Arrival  (ToA)  and  Received  Signal  Strength  (RSS)  based  mini¬ 
max  semidefinite  programming  (SDP)  algorithms  were  adopted  for  cooperative  localiza¬ 
tion  where  inter-target  communication  capability  was  exploited  for  coverage  extension  and 
accuracy  enhancement.  We  derived  the  Cramer-Rao  Lower  Bound  (CRLB)  to  verify  the 
benefits  of  cooperative  localization.  Various  semidefinite  algorithms,  including  Standard 
SDP  (SSDP),  Edge-based  SDP  (ESDP),  Node-based  SDP  (NSDP)  and  Component-wise  SDP 
(CSDP),  were  then  investigated  to  show  that  CSDP  was  the  best  among  them  in  terms  of 
localization  complexity  and  accuracy.  These  theoretical  analyses  had  all  been  corroborated 
with  simulations. 


H  From  Ranging  to  Localization:  Exploiting  Doppler  Effects 

Traditionally,  position  of  a  target  can  be  estimated  by  measuring  one  or  more  of  the  follow¬ 
ing  quantities:  angle-of-arrival  (AoA),  received  signal  strength  (RSS),  and  time-of-arrival 
(ToA).  For  AoA-based  techniques,  tracking  of  the  moving  target  relies  on  the  measurement 
of  the  angle  of  the  incoming  signal  arrival  at  each  anchor  node  (AN)  with  antenna  array. 
The  AoA  is  not  very  suitable  for  the  opportunistic  positioning,  because  the  antenna  array 
is  usually  large  and  expensive.  RSS-based  positioning  techniques  require  the  knowledge 
of  channel  path-loss  parameters  and  are  therefore  very  sensitive  to  the  estimation  of  these 
parameters.  Compared  to  AoA  and  RSS,  ToA  measurement  combined  with  the  triangula¬ 
tion  is  a  more  feasible  means  for  wireless  positioning  in  multipath  channel  environments, 
especially  when  broadband  signals  are  used  to  achieve  the  high  timing  resolution.  To  avoid 
the  need  of  the  common  time  base  between  the  target  and  each  anchor  node,  the  time  dif¬ 
ference  of  arrival  (TDoA)  technique  can  be  adopted.  However,  the  time  base  still  needs  to 
be  calibrated  among  all  anchor  nodes. 

In  order  to  avoid  these  problems  associated  with  the  traditional  methods,  we  propose 
a  new  localization  and  navigation  approach  based  on  the  Doppler  frequency  of  the  line- 
of-sight  (LoS)  path  between  the  moving  target  (MT)  and  each  AN.  As  long  as  the  AN  is 
able  to  identify  the  LoS  component  of  the  multipath  channel,  it  can  estimate  the  Doppler 
frequency  caused  by  the  movement  of  the  target.  After  the  Doppler  frequency  is  measured 
at  each  anchor  node,  the  radial  velocity  which  is  defined  as  the  relative  velocity  between 
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the  moving  target  and  each  anchor  node  can  be  estimated.  Suppose  that  we  have  N  ANs 
with  their  location  coordinates  known  to  the  system.  We  can  establish  N  equations  with 
the  radial  velocities  and  four  unknown  parameters:  the  moving  direction  of  the  target,  the 
absolute  velocity  of  the  target,  and  the  location  coordinates  of  the  moving  target.  Among 
them,  the  angle  of  moving  direction  and  the  absolute  velocity  are  two  nuisance  parameters. 
Solving  these  equations,  we  can  obtain  the  location  of  the  moving  target. 

Compared  to  existing  techniques,  the  hard¬ 
ware  cost  of  this  positioning  method  is  lower, 
because  no  antenna  array  is  needed.  In  addi¬ 
tion,  no  common  time  base  is  assumed  either  be¬ 
tween  the  moving  target  and  each  anchor  node 
or  among  the  anchor  nodes.  This  is  a  great  ad¬ 
vantage  in  comparison  with  the  conventional 
ToA  and  TDoA  techniques.  Besides,  our  tech¬ 
nique  benefits  from  the  motion  of  target,  perfor¬ 
mance  improves  as  velocity  increases. 

H.l  System  Description 

Suppose  we  have  N  anchor  nodes  (AN)  with  the 
location  coordinates  {{xi,yi)}^=l  of  the  ?'th  AN 
known  a  priori  to  the  system.  The  signal  transmission  can  be  either  from  MT  to  AN  or 
AN  to  MT.  Our  technique  can  be  applied  independent  of  the  direction  of  signal  flow.  The 
target  moves  in  a  direction  6  with  a  velocity  v  m/ sec.  Let  a*  be  the  angle  of  arrival  (AoA)  of 
the  signal  received  at  the  ith  AN,  as  shown  in  Fig.16. 

Then  the  radial  velocity  of  the  moving  target  with  respect  to  the  ?th  anchor  node  is 

Vi  =  V  COs(7T  +  CCj  —  6)  =  —v  cos(ccj  —  6)  (105) 

where  i  e  [1,  N],  ct*  G  [0,  2%),  and  0  e  [0,  2ir )  Our  objective  is  to  recover  the  location  coordi¬ 
nates  (xt,  yt)  of  the  moving  target  from  vt  and  { (x,.  yi)}f=x.  Solving  the  nonlinear  equation 
(105)  by  using  measurements  of  {vi}f=l,  we  can  estimate  the  location  of  the  moving  target. 
To  show  the  feasibility  of  this  approach,  we  first  present  a  very  direct  approach. 


Figure  16:  System  setup. 
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(a)  (b)  (c) 

Figure  17:  Error  of  mislocalization  (a)  1  anchor  node;  (b)  2  anchor  nodes;  (c)  3  anchor  nodes. 

H.1.1  Straightforward  Approach 

In  the  region  of  interest  where  MT  could  possibly  reside,  we  consider  every  possible  can¬ 
didate  location  (V,  y').  For  each  candidate  location,  calculate  the  radial  velocities  {u-}^  at 
all  ANs  assuming  velocity  v  and  direction  dare  known.  After  calculating  the  radial  veloci¬ 
ties  for  all  possible  locations,  we  choose  that  particular  location  which  gives  radial  velocities 
that  match  closest  to  the  actual  radial  velocities  {vi}^=1  obtained  from  the  doppler  spectrum. 
For  doing  this  we  define  the  error  of  mislocalization  given  by 

N 

E  =  I>-^)2 

2—1 

We  plot  the  error  of  mislocalization  as  shown  in  Fig.  17  for  N  =  (1,  2,  3)  anchor  nodes  that 
are  uniformly  located  on  the  circle  with  the  radius  of  2,  and  the  target  is  marked  by  the 
red  solid  triangle.  When  one  anchor  node  is  used,  the  target  can  be  on  either  of  two  half 
lines  where  the  error  is  0.  This  shows  the  approach  is  essentially  equivalent  to  the  AoA 
based  technique  without  any  antenna  array,  the  difference  being  that  here  we  have  two  half 
lines.  When  there  are  two  anchor  nodes,  the  target  can  be  found  at  the  intersection  of  the 
ambiguity  ranges  of  the  two  anchor  nodes.  Fig.  17  (b)  shows  that  the  error  surface  has  two 
minima  in  the  area  of  interest.  Only  one  corresponds  to  the  position  of  the  target.  The  other 
one  is  false.  As  shown  in  Fig.  17  (c),  when  three  anchor  nodes  are  used,  the  false  solution 
can  be  eliminated.  This  way,  the  unique  location  of  the  moving  target  can  be  estimated 
by  determining  the  minimum  peak  of  the  error  plot.  It  should  be  noted  that  {TyjO  ,  and 
(xt.  yt )  vary  as  the  target  moves.  However,  when  the  anchor  nodes  are  far  from  the  target. 
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these  parameters  can  be  assumed  to  be  constant  during  one  localization  processing  interval. 
Though  it  is  pretty  straightforward,  it  is  impractical  as  it  is  computationally  intensive  to  cal¬ 
culate  the  mislocalization  error  for  every  single  candidate  location  point  in  such  a  large  area 
covered  by  ANs.  Thus  we  should  explore  for  more  sophisticated  methods  for  practicality 


Our  goal  in  this  work  is  to  estimate  the  MT  location  co¬ 
ordinates  ( xt,iit )  by  solving  Eq.  (105).  As  for  the  solution  of 
this  nonlinear  equation,  no  linearization  technique  has  been 
found.  Since  among  all  the  unknown  parameters,  the  abso¬ 
lute  velocity  v  and  the  moving  direction  6  are  two  nuisance 
parameters,  if  at  least  one  or  both  of  v  and  6  are  known  by  the 
ANs,  the  problem  will  be  significantly  simplified.  Hence  in  the 
following  sections  we  present  the  solutions  to  the  localization 
equations  with  different  levels  of  prior  knowledge  of  v  and  6. 
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H.2  Case  I:  Known  v,  Known  9  / 

/ 

Consider  the  case  when  both  the  velocity  and  moving  direc¬ 
tion  of  the  target  are  known  at  the  ANs.  This  scenario  could  Figure  18:  Only  the  set  of  cor- 
occur  when  the  target  is  mounted  on  a  vehicle  that  is  moving  rect  AoA  intersect  at  the  MT 
with  a  constant  speed  in  the  same  direction  along  a  straight  location 
road.  The  velocity  of  the  target,  v  can  be  read  electronically 
from  the  speedometer  of  the  vehicle  and  the  moving  direction 
of  the  target,  6  can  be  obtained  from  the  gyroscope  on  the  vehicle. 


H.2.1  The  Basic  Algorithm 

Recollect  from  Eq.  (105),  when  both  v  and  6  are  known,  we  have 

cos  (a*  —  6)  =  — i  E  [1,  A].  (106) 

We  know  cc,  e  [0,  2tt)  and  9  e  [0,  27t),  and  hence  («;  —  9)  e  Lety*  =  cos_1(— ih/v)  be 

the  principal  solution.  Then  (a,  —  0)  is  one  of  the  four  possible  solutions  {y,,  2tt  —  y,,  -y,,  y,;  — 
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27r}  and  hence  a*  is  one  of  the  four  possible  solutions  {7*  +  0,  27r  —  7;  +  9 ,  —7 ,  +  9,9  —  2tt  +  7,;}. 

Using  the  known  condition  a*  e  [0,  2n),  we  can  directly  eliminate  two  among  four  pos¬ 
sible  solutions.  Now,  we  are  left  with  two  possible  solutions  of  the  AoA  a,.  Hence,  the 
problem  reduces  to  AoA  based  location  estimation.  Only  one  of  the  two  solutions  is  the 
correct  AoA.  The  correct  AoA  of  all  the  N  ANs  intersect  at  a  single  point  which  is  the  actual 
location  of  the  target.  The  other  angles  will  never  intersect  at  a  single  point  in  ideal  scenar¬ 
ios  where  the  noise  is  absent  (see  Fig.  18).  Thus  our  task  is  to  construct  half-lines  passing 
through  the  ith  AN  at  AoA  «,  and  select  the  particular  set  which  gives  single  intersect  point. 

H.2.2  Ambiguity  Issues 

Large  Scale  Ambiguity 

Ideally,  in  the  absence  of  noise,  the  correct  set  of  AoA  derived  at  each  anchor  node  intersects 
at  a  single  point,  the  actual  MT  location  (27,  yt ).  However  imperfect  estimation  in  the  pres¬ 
ence  of  noise  leads  to  deviated  AoAs  which  intersect  at  three  different  points  and  form  a 
triangle.  This  sometimes  results  in  a  false  set  of  AoAs  intersecting  to  form  a  smaller  triangle 
than  the  correct  set  of  AoAs  at  a  distance  very  far  away  from  actual  target  and  hence  the  lo¬ 
cation  error  shoots  up  suddenly  to  a  large  value  and  significantly  affect  the  performance  on 
average.  This  can  be  better  understood  from  the  error  mislocalization  function  with  perfect 
radial  velocities  replaced  by  the  estimated  vt  as  E  =  ~  vi )2- 

As  shown  in  Fig.  19(a)  for  three  anchor  nodes  located  at  1  km  apart,  multiple  minima 
occur.  It  can  be  observed  that  the  false  solution  occurs  very  far,  generally  lying  outside  the 
region  bounded  by  the  three  ANs.  Hence  in  order  to  avoid  the  false  solutions,  we  constrain 
our  region  of  interest  only  to  the  region  formed  by  the  set  of  ANs  assuming  that  only  those 
anchor  nodes  which  surround  the  MT  are  considered  for  estimation.  Another  way  of  elim¬ 
inating  the  false  minimum  is  to  consider  the  location  estimate  over  consecutive  intervals. 
The  false  minimum  will  change  drastically  while  the  correct  estimate  will  be  slowly  vary¬ 
ing.  We  need  to  pick  the  right  location  estimate  that  varies  slowly  thus  avoiding  the  false 
minimum.  This  problem  of  multiple  solutions  however  is  not  entirely  solved  by  constrain¬ 
ing  the  region  alone.  If  the  MT  and  any  pair  of  ANs  are  collinear,  then  the  multiple  solutions 
may  occur  even  for  a  constrained  region  as  shown  in  Fig.  19(b).  This  can  be  understood  in¬ 
tuitively  that  the  two  anchor  nodes  are  essentially  giving  the  same  information  about  the 
MT  location  and  not  giving  three  different  equations  to  solve.  To  tackle  this  problem,  we 
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(a)  (b) 

Figure  19:  (a)  Two  peaks  occur  in  error  plot  where  only  one  corresponds  to  actual  solution.  The 
false  solution  lies  outside  the  triangle  covered  by  the  three  anchor  nodes,  (b)  Two  peaks  occur  in 
error  plot  even  within  the  constrained  triangle  region 

introduced  at  least  one  more  anchor  nodes  which  introduces  more  equations.  Then  the  mul¬ 
tiple  solutions  will  merge  to  a  single  peak  thus  resolving  the  ambiguity. 


Small  Scale  Ambiguity 

In  the  presence  of  large  noise,  the  centroid  of  the  triangle  formed  by  joining  the  intersection 
points  of  AoA  half  lines  deviates  too  much  from  the  actual  location.  Sometimes,  the  centroid 
may  even  fall  outside  the  bounded  region.  As  a  remedy  we  introduce  AoA  sector  scanning 
where  we  scan  each  AoA  between  the  ranges  (AoA-Th)  to  (AoA+Th)  and  Th  is  the  possible 
error  that  occurs  in  AoA  as  shown  in  Fig.  20(a).  These  scanned  regions  of  all  ANs  will 
overlap  with  each  each  other.  The  MT  resides  in  the  region  where  the  AoA  sectors  of  all 
the  ANs  overlap.  We  then  choose  the  centroid  of  this  overlap  region  as  the  estimate  of  the 
location  of  MT.  This  technique  ensures  that  the  estimate  lies  within  the  bounded  region 
and  gives  an  estimate  closest  to  the  correct  solution.  To  obtain  the  AoA  error  threshold  Th, 
we  plot  the  histogram  of  AoA  error  as  shown  in  Fig.  20(b).  If  Th  is  too  large  then  we  are 
unnecessarily  considering  a  larger  overlapping  region  and  degrading  the  estimation  since 
for  most  of  the  cases  error  is  small.  On  the  other  hand,  if  Th  is  too  small  then  we  end  up 


57 


AN  3 


Histogram  of  error  in  AoA  at  SNR  =  30dB 


(a)  (b) 

Figure  20:  (a)  At  each  AN,  the  region  from  AoA-Th  to  AoA+Th  is  considered  and  the  MT  location 
falls  in  the  overlapping  region,  (b)  Histogram  of  Error  in  AoA  at  SNR  =  30dB.  The  threshold  is 
chosen  at  the  point  where  occurrence  of  error  falls  too  low. 

not  having  any  overlap  region  at  all  in  too  many  cases.  As  a  tradeoff  between  the  two 
extremes,  we  chose  the  threshold  Th  such  that  the  number  of  cases  with  AoA  error  greater 
than  Th  is  very  small  (say  <  2%)  as  most  of  the  error  is  concentrated  closer  to  smaller  error 
as  shown  in  Fig.  20(b).  When  the  error  is  greater  than  the  threshold,  we  average  over 
the  regions  overlapped  fewer  times  say  A  —  1,  A  —  2,  A  —  3  and  so  on  and  yet  performs 
better  than  smallest  triangle  method.  Please  note  that  the  threshold  varies  as  per  SNR  value. 
We  can  clearly  observe  in  Fig.  20(a)  that  this  method  gives  smaller  overlap  region  and 
closer  to  the  actual  location  than  the  triangle  formed  by  the  intersection.  Increasing  anchor 
nodes  further  reduces  this  region  of  overlap  and  thus  improving  the  performance.  Our 
AoA  sector  scanning  is  a  low  complexity  technique  to  determine  location  from  imperfect 
AoA  estimation.  There  are  several  other  existing  techniques  that  give  optimum  location 
estimate  (see  [9], [10], [11], [12])  that  can  be  applied  here. 

H.2.3  Algorithm  Flowchart 

Our  algorithm  to  estimate  the  location  of  MT  for  known  velocity  and  motion  direction  is 
summarized  in  the  flowchart  shown  in  Fig.  21. 
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Figure  21:  Flowchart  of  the  localization  technique  for  known  v  and  6  of  MT 

H.2.4  Simulations 

This  section  aims  to  analyze  the  performance  of  the  algorithm  and  study  the  influence  of  the 
number  of  anchor  nodes  and  temporal  averaging  on  the  performance.  Simulation  setup  in¬ 
volves  N  anchor  nodes  arranged  on  a  circle,  each  separated  from  its  neighbors  by  a  distance 
of  800m.  Depending  on  the  signal  strength  received  at  each  anchor  node,  we  can  further 
constrain  our  region  of  interest  thus  improving  accuracy.  As  an  example,  consider  the  setup 
as  shown  in  Fig.  22(a)  where  there  are  four  ANs.  If  the  received  signal  strength  is  greater  at 
one  anchor  node,  then  the  mobile  terminal  is  located  closer  to  that  anchor  node,  and  we  can 
constrain  our  region  of  interest  to  be  the  400m  x  400m  square  represented  by  shaded  region. 
We  are  right  now  considering  only  LOS  channel  model  and  the  velocity  components  are 
derived  from  the  doppler  shift  equations.  In  order  to  suppress  the  noise  influence  on  the 
performance,  the  measured  velocity  component  is  averaged  M  times  temporally  during  ev¬ 
ery  localization  interval.  With  this  setup  and  estimation  method,  the  simulations  have  been 
done  and  the  root  mean  square  error  (RMSE)  plots  have  been  obtained  for  velocity  range 
of  mobile  terminal  36kmph-90kmph  (10mps-25mps)  for  N  =  4,8  and  M  =  128  as  shown  in 
Fig.  22(b).  We  can  clearly  see  that  the  performance  improves  as  the  number  of  anchor  nodes 
are  increased. 
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AN  1  AN  1 


(a)  (b) 

Figure  22:  (a)  Setup  of  4  anchor  nodes  separated  by  800m.  MT  is  closer  to  AN  1.  Hence  the  region 
of  interest  is  constrained  to  the  400m  x  400m  shaded  area,  (b)  RMSE  vs  SNR  plot  for  known  direction 
and  velocity.  Performance  improves  as  number  of  anchor  nodes  increased  and  when  averaged  over 
more  intervals  of  time 

H.3  Case  II:  Known  v,  Unknown  9 

For  this  case,  the  scenario  is  very  similar  to  the  first  one  except  that  the  moving  direction 
information  of  the  target  does  not  need  to  be  measured.  Since  the  direction  9  is  unknown, 
the  method  in  the  preceding  section  can  not  be  applied  here.  Thus  we  should  aim  to  elim¬ 
inate  the  nuisance  parameter  9  and  solve  the  nonlinear  equations  by  relating  the  required 
quantities  location  coordinates  (xt,  yt)  with  the  known  quantities  v,  vt .  and  {xt, 

H.3.1  The  Basic  Algorithm 

To  derive  the  relationship  between  the  MT  location  co-ordinates  (xt,  yt)  and  the  AoA  ay,  at 
the  ith  anchor  node  we  consider  the  relation  tan(ay)  =  Let 

A  =  tan”1  ( (107) 

\  X  -  Xi ) 

We  can  easily  derive  the  expression  for  ay  as 

f  fa,  if  ay  <  |  ; 

ai  =  l  t t  +  Pi,  if  |  <  ay  <  ^  ;  (108) 

[  27T  +  Pi,  if  ^  <  ai  <  2vr 
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Using  the  identity  tan(/c7r  +  ui)  —  tan(ca),  we  get 


tan  —  9)  =  tan  (/%  —  6*) 


(109) 


To  link  known  quantities  v,  vtJ  and  {x,  ,  yt  ]U ,  with  the  AoA  a*,  we  can  re-express  equation 
(106)  as 


1  +  tan2(oj  —  9)  = 


which  implies  that 


tan  (a,  —  9)  =  —  —  1 


Let  Tji  =  tan  ,  /  h  -  1  be  the  principal  solution.  Then, 


tan(oj  —  6)  =  ±tan(r/j) 


(110) 


Equating  (109)  and  (110),  we  get 


tan(/?j  —  9)  =  ±  tan(77j) 


(111) 


In  order  to  eliminate  0,  consider  the  AN  pair  ( i,j )  and  perform  the  subtraction  to  get 


tan (fli  —  fij)  =  tan(±r7i  ±  rjj). 


(112) 


In  (112),  the  right  hand  side  (RHS)  can  be  calculated  from  v  and  vl  while  for  the  left  hand 
side  (LHS),  using  the  definition  of  /?*  given  in  (107)  and  the  known  identity,  tan  (A  —  B)  = 


tan  A— tan  B 
l+tanAtanB 


,  we  get 


ta n(/3i-/3j)=  (a  _  +  tt») 

v  V  X—Xi  X — Xj  '  '  V  X—XiX—Xj' 


(113) 


Substituting  LHS  and  RHS  back  in  (112)  and  reordering  the  terms,  we  get  the  equation 


(*  -  Aid)2  +  (v~  Bij?  =  Al  •  +  Bfj  -  Cij 


(114) 


where 


a  _  X%  T  Xj  l)i  Xjj 

2  2  tan  (±77*  ±  rjj) 


d  y*  +  Vi 
B.,i  =  —5— 


<2/  j  oc  ^ 


2  tan(±-//j  ±  rjj) 


Cij  =  XiXj  +  //,//,  + 


XiDj  %j]Ji 

tan(±77j  ±  rjj) 
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The  RHS  in  (114)  on  simplification  gives 


Ah  +  Bh  -  C,<  =  A  I  1 


1,3  '  1,3  1,3  4  V  '  tan2(±77i  ±  rjj) 


>  0 


(115) 


where  is  the  distance  between  the  AN  i  and  AN  j.  Thus  (114)  is  of  the  form  of  equation 
of  a  circle  with  center  co-ordinates  as  (Ajj,  Bhj)  and  radius  as 


r  “  (/A  +  Bh  -  C‘.i  -  +  tan2(±I)i  ±Vjy 


Interestingly,  these  circles  have  the  following  properties: 


•  The  circle  corresponding  to  the  AN  pair  (i,  j),  ( x  —  Aitj)2  +  (y  —  Bij)'2  =  A2  -  +  B2  -  —  ChJ 
passes  through  ANs  i  and  j.  This  can  be  verified  easily  by  substituting  the  AN's  co¬ 
ordinates  in  the  equation. 


•  For  N  ANs  we  get  (A  —  1)  independent  circles  for  each  of  2N  different  combinations 
of  ±77 ,  of  which  only  one  particular  combination  of  (A  —  1)  circles  intersect  at  a  single 
point  (other  than  at  AN  itself)  which  is  the  MT  location. 


H.3.2  Ambiguity  Issues 

Similar  to  the  known  velocity  (v)  and  known  direction  (0)  case,  ambiguity  issues  arise  in 
this  unknown  direction  case  due  to  imperfect  estimation  of  radial  velocities  at  the  ANs. 
The  large  scale  ambiguity  is  combated  by  constraining  the  region  to  the  area  bounded  by 
ANs  and  by  increasing  the  number  of  ANs.  The  small  scale  ambiguity  is  combatted  in 
similar  way  as  AoA  sector  scanning,  except  that  instead  of  AoA  we  have  the  quantity  Ptj  = 
tan(±r/,  ±  r/j )  and  we  plot  error  histogram  from  Pt  J  and  sweep  it  between  (Pl  ;]  ±  T h)  for  each 
circle  and  determine  the  overlap  region.  The  centroid  of  this  overlapping  region  is  taken  as 
the  location  estimate. 


H.3.3  Algorithm  Flowchart 

Our  algorithm  to  estimate  location  of  MT  for  known  velocity  and  unknown  direction  is 
summarized  in  the  form  of  flowchart  shown  in  Fig.  23. 

H.3.4  Simulations 
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Figure  23:  Flowchart  of  the  localization  technique  for  known  v  and  unknown  6  of  MT 


The  simulation  setup  is  exactly  same  as  the  known  di¬ 
rection  case  as  explained  in  detail  earlier.  The  perfor¬ 
mance  plot  is  as  shown  in  the  Fig.  24.  The  performance 
of  our  technique  for  the  unknown  direction  case  is  sim¬ 
ilar  to  that  of  the  known  direction  case  but  this  is  ob¬ 
tained  at  the  cost  of  increased  computation  and  mini¬ 
mum  number  of  ANs  required. 


H.4  Summary 


Figure  24:  RMSE  vs  SNR  plot  for  the 
unknown  direction  and  known  veloc¬ 
ity.  Performance  improves  as  number 

of  anchor  nodes  increased.  The  per- 
We  presented  a  new  approach  to  find  the  location  of  a  formance  for  known  and  unkn0wn 

moving  target  using  the  observed  radial  velocities  from  djrecti0n  cases  is  similar 
Doppler.  We  specifically  focused  on  two  cases  -  1)  Ve¬ 
locity  ( v )  and  direction  of  motion  (0)  of  the  moving  tar¬ 
get  are  known  2)  Only  velocity  ( v )  is  known  and  direc¬ 
tion  (6)  is  unknown.  We  provided  methods  to  estimate  location  coordinates  of  the  moving 
target  by  solving  a  set  of  nonlinear  equations.  Our  approach  imposes  no  requirement  of 
common  time  base  between  any  two  nodes  and  avoids  large  expensive  antenna  arrays  at 
the  anchor  nodes.  Since  our  technique  benefits  from  the  motion  of  the  target,  the  perfor¬ 
mance  improves  as  the  mobility  increases. 
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